Signal Classification under structure sparsity constraints

TL;DR

This paper develops novel sparse signal processing algorithms for object classification, demonstrating their effectiveness across medical imaging, radar detection, and general image categorization tasks.

Contribution

It introduces new optimization formulations for structured sparse representations and provides tractable solutions, applying them to diverse real-world classification problems.

Findings

Effective classification in medical imaging

Successful detection of buried objects with radar

Improved accuracy in general object recognition

Abstract

Object Classification is a key direction of research in signal and image processing, computer vision and artificial intelligence. The goal is to come up with algorithms that automatically analyze images and put them in predefined categories. This dissertation focuses on the theory and application of sparse signal processing and learning algorithms for image processing and computer vision, especially object classification problems. A key emphasis of this work is to formulate novel optimization problems for learning dictionary and structured sparse representations. Tractable solutions are proposed subsequently for the corresponding optimization problems. An important goal of this dissertation is to demonstrate the wide applications of these algorithmic tools for real-world applications. To that end, we explored important problems in the areas of: 1. Medical imaging: histopathological images acquired from mammalian tissues, human breast tissues, and human brain tissues. 2. Low-frequency (UHF to L-band) ultra-wideband (UWB) synthetic aperture radar: detecting bombs and mines buried under rough surfaces. 3. General object classification: face, flowers, objects, dogs, indoor scenes, etc.

Figures

Tables

Table 1. Table 1 : Complexity analysis for different dictionary learning methods.

Method	Complexity	Running time
DFDL	$c^{2}kN\left(2d+L^2\right)$	$\sim$ 0.5 hours
LC-KSVD[37]	$c^{2}kN\left(2d+2ck+L^2\right)$	$\sim$ 3 hours
Nayak’s et al.[48](∗)	$c^{2}kN\left(2d+2qck\right)+c^{2}d{k}^2$	$\sim$ 8 hours
FDDL[38](∗)	$c^{2}kN\left(2d+2qck\right)+c^{3}d{k}^2$	$>$ 40 hours

Table 2. Table 2 : Estimated number of operations required in different dictionary learning methods.

Method	$q=1$	$q=3$	$q=10$
DFDL	$6.6\times {10}^{10}$	$6.6\times {10}^{10}$	$6.6\times {10}^{10}$
LC-KSVD[37]	$1.06\times {10}^{11}$	$1.06\times {10}^{11}$	$1.06\times {10}^{11}$
Nayak’s et al.[48]	$8.92\times {10}^{10}$	$1.692\times {10}^{11}$	$4.492\times {10}^{11}$
FDDL[38]	$9.04\times {10}^{10}$	$1.704\times {10}^{11}$	$4.504\times {10}^{11}$

Table 3. Table 3 : CONFUSION MATRIX: IBL.

Class	UDH	DCIS	Method
UDH	91.75	8.25	WND-CHARM(∗) [53]
	68.00	32.00	SRC(∗) [1]
	93.33	6.67	SHIRC [20]
	84.80	15.20	FDDL [38]
	90.29	9.71	LC-KSVD [37]
	85.71	14.29	Nayak’s et al.[48]
	96.00	4.00	DFDL
DCIS	5.77	94.23	WND-CHARM(∗) [53]
	44.00	56.00	SRC(∗) [1]
	10.00	90.00	SHIRC [20]
	10.00	90.00	FDDL [38]
	14.86	85.14	LC-KSVD [37]
	23.43	76.57	Nayak’s et al.[48]
	0.50	99.50	DFDL

Table 4. Table 4 : CONFUSION MATRIX: ADL ( % percent \% ).

	Kidney		Lung		Spleen
Class	H.	I.	H.	I.	H.	I.	Method
H.	83.27	16.73	83.20	16.80	87.23	12.77	WND-CHARM(∗) [53]
	87.50	12.50	72.50	27.50	70.83	29.17	SRC(∗) [1]
	82.50	17.50	75.00	25.00	65.00	35.00	SHIRC [20]
	83.26	16.74	93.15	6.85	86.94	13.06	FDDL [38]
	86.84	13.16	85.59	15.41	89.75	10.25	LC-KSVD [37]
	73.08	26.92	89.55	10.45	86.44	13.56	Nayak’s et al.[48]
	88.21	11.79	96.52	3.48	92.88	7.12	DFDL
I.	14.22	85.78	14.31	83.69	10.48	89.52	WND-CHARM(∗) [53]
	25.00	75.00	24.17	75.83	20.83	79.17	SRC(∗) [1]
	16.67	83.33	15.00	85.00	11.67	88.33	SHIRC [20]
	19.88	80.12	10.00	90.00	8.57	91.43	FDDL [38]
	19.25	81.75	10.89	89.11	8.57	91.43	LC-KSVD [37]
	26.92	73.08	25.90	74.10	6.05	93.95	Nayak’s et al.[48]
	9.92	90.02	2.57	97.43	7.89	92.01	DFDL

Table 5. Table 5 : CONFUSION MATRIX: TCGA ( % percent \% ).

Class	Not MVP	MVP	Method
Not VMP	76.68	23.32	WND-CHARM[53]
	92.92	7.08	Nayak’s et al.[48]
	96.46	3.54	LC-KSVD[37]
	92.04	7.96	FDDL[38]
	94.69	5.31	DFDL
MVP	21.62	78.38	WND-CHARM[53]
	16.22	83.78	Nayak’s et al.[48]
	8.10	91.90	LC-KSVD[37]
	18.92	81.08	FDDL[38]
	5.41	94.59	DFDL

Table 6. Table 6 : Complexity analysis for proposed efficient algorithms and their original versions

Method	Complexity	Plugging numbers
O-DLSI-D	$Cqk{d}^3$	$6.25\times {10}^{12}$
E-DLSI-D	$C{d}^3+Cqdk\left(qk+k\right)$	$2.52\times {10}^{10}$
O-FDDL-X	$C^{2}k\left(dn+qCkn+Cdk\right)$	$1.51\times {10}^{11}$
E-FDDL-X	$C^{2}k\left(dn+qCnk+dk\right)$	$1.01\times {10}^{11}$
O-FDDL-D	$Cdk\left(qk+C^{2}n\right)$	${10}^{11}$
E-FDDL-D	$Cdk\left(Cn+Cqk\right)+C^3{k}^{2}n$	$2.8\times {10}^{10}$

Table 7. Table 7 : Complexity analysis for different dictionary learning methods

Method	Complexity	Plugging numbers
O-DLSI	$Ck\left(kd+dn+qkn\right)+Cqk{d}^3$	$6.25\times {10}^{12}$
E-DLSI	$Ck\left(kd+dn+qkn\right)+$ $C{d}^3+Cqdk\left(qk+d\right)$	$3.75\times {10}^{10}$
O-FDDL	$C^{2}dk\left(n+Ck+Cn\right)+$ $+C{k}^{2}q\left(d+C^{2}n\right)$	$2.51\times {10}^{11}$
E-FDDL	$C^{2}k\left(\left(q+1\right)k\left(d+Cn\right)+2dn\right)$	$1.29\times {10}^{11}$
O-COPAR	$C^3{k}^2\left(2d+Ck+qn\right)+Cqk{d}^3$	$6.55\times {10}^{12}$
E-COPAR	$C^3{k}^2\left(2d+Ck+qn\right)+$ $+C{d}^3+Cqdk\left(qk+d\right)$	$3.38\times {10}^{11}$
LRSDL	$C^{2}k\left(\left(q+1\right)k\left(d+Cn\right)+2dn\right)$ $C^{2}dkn+\left(q+q_2\right)d{k}^2$	$1.3\times {10}^{11}$

Table 8. Table 8 : Overall accuracy (mean ± plus-or-minus \pm\leavevmode\nobreak\ standard deviation) (%) of different dictionary learning methods on different datasets. Numbers in parentheses are number of training samples per class.

	Ext. YaleB (30)	AR (20)	AR gender (250)	Oxford Flower (60)	Caltech 101 (30)	COIL100 (10)
SRC [1]	97.96 $\pm$0.22	97.33 $\pm$0.39	92.57 $\pm$0.00	75.79 $\pm$0.23	72.15 $\pm$0.36	81.45 $\pm$ 0.80
LC-KSVD1 [37]	97.09 $\pm$0.52	97.78 $\pm$0.36	88.42 $\pm$1.02	91.47 $\pm$1.04	73.40 $\pm$0.64	81.37 $\pm$0.31
LC-KSVD2 [37]	97.80 $\pm$0.37	97.70 $\pm$0.23	90.14 $\pm$0.45	92.00 $\pm$0.73	73.60 $\pm$0.53	81.42 $\pm$0.33
DLSI [43]	96.50 $\pm$0.85	96.67 $\pm$1.02	93.86 $\pm$0.27	85.29 $\pm$1.12	70.67 $\pm$0.73	80.67 $\pm$0.46
DLRD [79]	93.56 $\pm$1.25	97.83 $\pm$0.80	92.71 $\pm$0.43	-	-	-
FDDL [4]	97.52 $\pm$0.63	96.16 $\pm$1.16	93.70 $\pm$0.24	91.17 $\pm$0.89	72.94 $\pm$0.26	77.45 $\pm$1.04
$D^2{L}^2{R}^2$[45]	96.70 $\pm$0.57	95.33 $\pm$1.03	93.71 $\pm$0.87	83.23 $\pm$1.34	75.26 $\pm$0.72	76.27 $\pm$0.98
COPAR [3]	98.19 $\pm$0.21	98.50 $\pm$0.53	95.14 $\pm$0.52	85.29 $\pm$0.74	76.05 $\pm$0.72	80.46 $\pm$0.61
JDL [33]	94.99 $\pm$0.53	96.00 $\pm$0.96	93.86 $\pm$0.43	80.29 $\pm$0.26	75.90 $\pm$0.70	80.77 $\pm$0.85
JDL∗ [33]	97.73 $\pm$0.66	98.80 $\pm$0.34	92.83 $\pm$0.12	80.29 $\pm$0.26	73.47 $\pm$0.67	80.30 $\pm$1.10
SRRS [80]	97.75 $\pm$0.58	96.70 $\pm$1.26	91.28 $\pm$0.15	88.52 $\pm$0.64	65.22 + 0.34	85.04 $\pm$0.45
LRSDL	98.76 $\pm$0.23	98.87 $\pm$0.43	95.42 $\pm$0.48	92.58 $\pm$0.62	76.70 $\pm$0.42	84.35 $\pm$0.37

Table 9. Table 9 : Training and test time per sample (seconds) of different dictionary learning method on the Oxford Flower dataset ( n c = 60 , d = 10000 , C = 17 , K ≈ 40 × 17 formulae-sequence subscript 𝑛 𝑐 60 formulae-sequence 𝑑 10000 formulae-sequence 𝐶 17 𝐾 40 17 n_{c}=60,d=10000,C=17,K\approx 40\times 17 ).

	SRC	LCKSVD1	LCKSVD2	DLSI	FDDL	$D^2{L}^2{R}^2$	COPAR	JDL	SRRS	LRSDL
Train	0	1.53e3	1.46e3	1.4e3	9.2e2	$>$1 day	1.8e4	7.5e1	3.2e3	1.8e3
Test	3.2e-2	6.8e-3	6.1e-3	4e-3	6.4e-3	3.3e-2	5.5e-3	3.6e-3	3.7e-3	2.4e-2

Table 10. Table 10 : Overall average classification (%) of different methods on different polarization combinations, with or without the non-negativity constraint (NN).

	Separate targets					All targets
	VV	HH	VV+HH	HH+HV	ALL	VV	HH	VV+HH	HH+HV	ALL
SVM	62.17	74.47	65.00	77.77	70.33	83.68	89.36	83.36	66.36	79.78
CR	76.32	83.43	86.26	85.93	82.95	86.22	91.64	87.14	91.28	90.64	wo/ NN
CC	76.32	83.43	79.47	86.21	79.93	86.22	91.64	87.96	92.98	81.80
SM	76.32	83.43	81.07	86.52	85.87	86.22	91.64	89.66	95.12	94.14
GT	69.79	78.79	77.83	85.12	81.75	81.90	88.40	87.28	91.86	90.46
CR	79.05	89.41	79.85	88.52	84.92	86.32	92.98	87.48	94.86	94.26	w/ NN
CC	79.05	89.41	83.99	92.27	85.11	86.32	92.98	89.20	93.82	87.08
SM	79.05	89.41	85.62	90.85	89.36	86.32	92.98	90.00	96.82	96.00
GT	75.55	88.97	81.85	90.57	87.28	84.16	90.32	88.10	94.58	94.12

Table 11. Table 11 : Comparison of different discriminative sparsity frameworks (more ticks represent better performance).

	SRC[1]	DLSI[43]	FDDL[4]	COPAR[3]	LRSDL
training speed	no	✓	-	✓	✓✓
shared features	-	-	-	✓	✓
low-rank constraint	-	-	-	-	✓
structured coefficients	-	-	✓	✓	✓
structured dictionaries	-	✓	✓	✓	✓
classification performance	-	-	✓	✓	✓✓

Equations

$$\hat{\mathbf{x}}=\arg \underset{\mathbf{x}}{\min }\{\Vert \mathbf{y}-\mathbf{Dx}\Vert +\lambda \Vert \mathbf{x}{\Vert }_1\}$$

$$\text{identity}(\mathbf{y})=\arg \underset{i}{\min }\{e_i\}$$

$$(\mathbf{D}^{*},\mathbf{X}^{*})=\arg\min_{\mathbf{D},\mathbf{X}}\frac{1}{N}\sum_{i=1}^{N}\Big{\{}\left\|\mathbf{y}_{i}-\mathbf{D}\mathbf{x}_{i}\right\|_{2}^{2}+\lambda\left\|\mathbf{x}_{i}\right\|_{1}\Big{\}}$$

$$(\mathbf{D}^{*},\mathbf{X}^{*})=\arg\min_{\mathbf{D},\mathbf{X}}\Big{\{}\left\|\mathbf{Y}-\mathbf{D}\mathbf{X}\right\|_{F}^{2}+\lambda\left\|\mathbf{X}\right\|_{1}\Big{\}}$$

$${\mathbf{X}}^{\ast }=\arg \underset{\mathbf{X}}{\min }\{\Vert \mathbf{Y}-\mathbf{DX}{\Vert }_F^2+\lambda \Vert \mathbf{X}{\Vert }_1\}$$

$${\mathbf{D}}^{\ast }$$

$$=$$

$${\mathbf{u}}_j$$

$${\mathbf{d}}_j$$

$$(\hat{\mathbf{D}},\hat{\mathbf{X}},\hat{\mathbf{A}})$$

$$\Vert {\mathbf{x}}_i{\Vert }_0\le L$$

$$\mathbf{x}=\arg \underset{\mathbf{x}}{\min }\Vert \mathbf{y}-\hat{\mathbf{D}}\mathbf{x}{\Vert }_2^2\text{leavevmode}\text{leavevmode}\text{subject to }\text{leavevmode}\Vert \mathbf{x}{\Vert }_0\le L$$

$$\text{identity}(\mathbf{y})=\arg \underset{j}{\max }(l=\hat{\mathbf{W}}\mathbf{x})$$

$$J_{\mathbf{D},\mathbf{X}}=\arg\min_{\mathbf{D},\mathbf{X}}\Big{\{}\sum_{i=1}^{c}r(\mathbf{Y}_{i},\mathbf{D},\mathbf{X}_{i})+\lambda_{1}\|\mathbf{X}\|_{1}+\lambda_{2}\big{(}\textrm{trace}(S_{W}(\mathbf{X})-S_{B}(\mathbf{X}))+\eta\|\mathbf{X}\|_{F}^{2}\big{)}\Big{\}}$$

$$\underset{\Vert {\mathbf{x}}_l{\Vert }_0\le L_i}{\min }\Vert {\mathbf{y}}_l-{\mathbf{D}}_i{\mathbf{x}}_l{\Vert }_2^2,\text{leavevmode}\text{leavevmode}\forall l=1,2,\dots ,N_i$$

$$\underset{\Vert {\bar{\mathbf{x}}}_m{\Vert }_0\le L_i}{\min }\Vert {\bar{\mathbf{y}}}_m-{\mathbf{D}}_i{\bar{\mathbf{x}}}_m{\Vert }_2^2,\text{leavevmode}\text{leavevmode}\forall m=1,2,\dots ,{\bar{N}}_i$$

$$\text{intra-class differences: }\frac{1}{N_i}\underset{\Vert {\mathbf{X}}_i{\Vert }_0\le L_i}{\min }\Vert {\mathbf{Y}}_i-{\mathbf{D}}_i{\mathbf{X}}_i{\Vert }_F^2$$

$$\text{inter-class differences: }\frac{1}{\overline{N_i}}\underset{\Vert to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}\mathbf{Y}}_i{\Vert }_0\le L_i}{\min }\Vert to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}\mathbf{X}}_i-{\mathbf{D}}_{i}to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}\mathbf{X}}_i{\Vert }_F^2$$

$$\mathbf{D}^{*}=\arg\min_{\mathbf{D}}\Big{(}\frac{1}{N}\min_{\|\mathbf{X}\|_{0}\leq L}\|\mathbf{Y}-\mathbf{D}\mathbf{X}\|_{F}^{2}-\frac{\rho}{\bar{N}}\min_{\|\hbox to0.0pt{\hskip 0.12498pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathbf{X}$}}\|_{0}\leq L}\|\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathbf{Y}$}}-\mathbf{D}\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color$\overline{\hbox{}}$}\hss}{\leavevmode\hbox{\set@color$\mathbf{X}$}}\|_{F}^{2}\Big{)},$$

$$({\mathbf{D}}^0,{\mathbf{X}}^0)=\arg \underset{\mathbf{D},\mathbf{X}}{\min }\{\Vert \mathbf{Y}-\mathbf{DX}{\Vert }_F^2+\lambda \Vert \mathbf{X}{\Vert }_1\},$$

$$L\approx \frac{1}{N}\sum _{i=1}^N\Vert {\mathbf{x}}_i^0{\Vert }_0.$$

$$\hat{\mathbf{x}}=\arg\min_{\mathbf{x}}\big{\{}\|\mathbf{y}-\mathbf{D}_{total}\mathbf{x}\|_{2}^{2}+\gamma\|\mathbf{x}\|_{1}\big{\}},$$

$$r_i(\mathbf{y})=\Vert \mathbf{y}-{\mathbf{D}}_i{\delta }_i(\hat{\mathbf{x}}){\Vert }_2$$

$${\mathbf{X}}^{\ast }=\arg \underset{\Vert \mathbf{X}{\Vert }_0\le L}{\min }\Vert \mathbf{Y}-\mathbf{DX}{\Vert }_F^2;\text{leavevmode}\text{leavevmode}to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}{\text{leavevmode}\text{set@color}\mathbf{X}}^{\ast }=\arg \underset{\Vert to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}\mathbf{X}{\Vert }_0\le L}{\min }\Vert to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}\mathbf{Y}-\mathbf{D}to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}\mathbf{X}{\Vert }_F^2.$$

$${\hat{\mathbf{X}}}^{\ast }=\arg \underset{\Vert \hat{\mathbf{X}}{\Vert }_0\le L}{\min }{\Vert \hat{\mathbf{Y}}-\mathbf{D}\hat{\mathbf{X}}\Vert }_F^2.$$

$${\mathbf{D}}^{\ast }$$

$$\displaystyle=\arg\min_{\mathbf{D}}\big{\{}-2\textrm{trace}(\mathbf{E}\mathbf{D}^{T})+\textrm{trace}(\mathbf{D}\mathbf{F}\mathbf{D}^{T})\big{\}}.$$

$$\mathbf{E}=\frac{1}{N}\mathbf{Y}{\mathbf{X}}^T-\frac{\rho }{\bar{N}}\bar{\mathbf{X}}{\bar{\mathbf{X}}}^T;\mathbf{F}=\frac{1}{N}\mathbf{X}{\mathbf{X}}^T-\frac{\rho }{\bar{N}}\bar{\mathbf{X}}{\bar{\mathbf{X}}}^T.$$

$${\mathbf{D}}^{\ast }=\arg \underset{\mathbf{D}}{\min }\left(\frac{1}{N}\underset{\Vert \mathbf{X}{\Vert }_0\le L}{\min }\Vert \mathbf{Y}-\mathbf{DX}{\Vert }_F^2-\frac{\rho }{\bar{N}}\underset{\Vert to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}\mathbf{X}{\Vert }_0\le L}{\min }\Vert to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}\mathbf{Y}-\mathbf{D}to0.0pt\text{leavevmode}\text{set@color}\overline{}\text{hss}\text{leavevmode}\text{set@color}\mathbf{X}{\Vert }_F^2\right),$$

$$\mathbf{E}=\frac{1}{N}\mathbf{Y}{\mathbf{X}}^T-\frac{\rho }{\bar{N}}\bar{\mathbf{X}}{\bar{\mathbf{X}}}^T;\mathbf{F}=\frac{1}{N}\mathbf{X}{\mathbf{X}}^T-\frac{\rho }{\bar{N}}\bar{\mathbf{X}}{\bar{\mathbf{X}}}^T.$$

$$\mathbf{D}^{*}=\arg\min_{\mathbf{D}}\Big{\{}-2\textrm{trace}(\mathbf{E}\mathbf{D}^{T})+\textrm{trace}\Big{(}\mathbf{D}\big{(}\mathbf{F}-\lambda_{\min}(\mathbf{F})\mathbf{I}\big{)}\mathbf{D}^{T}\Big{)}\Big{\}}$$

$$\text{subject to:}\Vert {\mathbf{d}}_i{\Vert }_2^2=1,i=1,2,\dots ,k.$$

$$\mathbf{D}^{*}=\arg\min_{\mathbf{D}}\{-2\textrm{trace}(\mathbf{E}\mathbf{D}^{T})+\textrm{trace}\big{(}\mathbf{D}(\mathbf{F}-\lambda_{\min}(\mathbf{F})\mathbf{I}_{k})\mathbf{D}^{T}\big{)}\}$$

Full text

Signal classification under

structured sparsity constraints

Tiep Huu Vu

\ChTitleVar\dept

Electrical Engineering \degreedateMay 2019 \honorsdegreeinfofor a baccalaureate degree

in Engineering Science

with honors in Engineering Science \documenttypeDissertation \submittedtoThe Graduate School \numberofreaders4 \honorsadviserHonors P. Adviser \secondthesissupervisorSecond T. Supervisor \honorsdeptheadDepartment Q. Head \advisor[Dissertation Advisor, Chair of Committee] Vishal Monga Associate Professor of Electrical Engineering \readerone[] William E. Higgins Distinguished Professor of Electrical Engineering \readertwo[] Kenneth Jenkins Professor of Electrical Engineering \readerthree[] Robert T. Collins Associate Professor of Computer Science and Engineering \readerfour[] Kultegin Aydin Professor of Electrical Engineering and Department Head

\psutitlepage

\psucommitteepage

\thesisabstract

Abstract

\thesistableofcontents

\thesislistoffigures

\thesislistoftables

\thesisacknowledgments

Acknowledgments

\thesismainmatter

Introduction

Motivation

Object classification is one of the most important problems in machine intelligence systems. This is a well-studied problem in many important domains of data processing today. Typically, we are given images belonging to two or more classes, and the challenge is to develop an effective procedure of determining identity of a new sample. Application of object classification ranges broadly from spam email detection, fraud transaction detection to fingerprint verification, cancer detection, face identification, military vehicle recognition, self-driving cars, etc (see Fig. 1).

In general, a classification system comprises of two main components: a feature extraction tool and a classifier. The feature extraction part plays a role in generating discriminative characteristics of each class. Traditional techniques in extracting features are hand-designed, such as Haar wavelets, Difference of Gaussians (DOG) filter, Gabor filters, histogram of oriented gradient (HOG) descriptors, SIFT descriptors[6], spatial pyramid matching[7] and many others. In some particular cases, the feature extraction also includes data dimension reduction step, which reduces the complexity of the classifier and hence, maximizes speed of the system. The classifier, which is often a result of a learning scheme, takes the extracted features as inputs and predicts identity of signals. Several classifier design techniques have been applied in practical applications, such as Naive Bayes [8], Decision Trees [9], Logistic Regression [10], Neural Networks [11], Support Vector Machines [12], and Deep Learning [13, 14] recently.

One classifier can be seen as a function $f$ that takes an input signal $\mathbf{x}$ (discrete or continuous) and generates an discrete output $y\in\mathcal{C}$ representing category/class of $\mathbf{x}$. Based on several training pairs $\{(\mathbf{x}_{1},y_{1}),(\mathbf{x}_{2},y_{2}),\dots,(\mathbf{x}_{N},y_{N})\}$, a classification algorithm tries to ‘learn’ the mapping $f$ such that $y_{n}\approx f(\mathbf{x}_{n})$ for as many as possible $n\in\{1,2,\dots,N\}$. A good classifier is one that not only could well capture relationship of training pairs, i.e. what it has seen, but also need to well represent new samples. This property, called generalization, is one of the most important factors in a classification system. However, attaining this property is usually a challenging task, as real-world classification problems exhibit difficulties. In particular, this dissertation emphasizes on addressing the following challenges:

• •

Many practical problems encounter the issue of insufficient training. From probability perspective, the identity of a signal can be formulated as a maximum likelihood estimation problem, i.e. estimating $p(y|\mathbf{x})$ for each possible category $y\in\mathcal{C}$. Class decisions are made using empirical estimates of the true densities learned from available training data. Without prior knowledge of the (limited) training data, a system can severely misestimate the density, resulting in a poor classifier with lack of generalization property. This problem becomes more severe when data dimensionality increases. In this case, the volume of the space drastically increases. In order to achieve a stable system, the amount of data required to support the result often grows exponentially with the dimensionality. Unfortunately, this high-dimensional data phenomenon is commonly encountered in signal classification problems. For instance, a small gray image of size $200\times 200$ has dimension of 40000. Particular images, such as hyperspectral or medical, have even more dimensions with million of pixels per channel and multiple channels per image. The problem persists even when low-dimensional image features are considered instead of entire images. Training insufficiency and high dimensionality are thus important concerns in signal classification.

• •

Another challenge is that in practical problems, generous training is often not only limited but also taken under presence of noise or occlusion. This might be caused by limitations of devices, e.g. optical sensors, microscopes, or radar transmitter/receiver, that capture signals. Some unexpected sources of noise or distortions are often incorporated into the classification frameworks. In a face identification system, training faces might be captured under a good light conditions and/or well-aligned; a test face image, however, is usually captured under bad conditions with low light, occlusion by other objects, or at different angles. In a radar system, captured signals are commonly heavily interfered by other signals in the field.

• •

Additionally, classification systems often face the problem of high within-class variance and small between-class variance. In a face identity problem, an image could include the face of the same individual wearing sunglasses, with hat, or at different angles. In a flower classification problem, a flower may be captured with different background, at different times of a day, or at different stages of the blooming process. This variability of samples from the same class is also called high within-class variance. Another challenge could happen in a classification problem is that different classes could share common patterns. Objects of different classes might be captured under the same background conditions. In histopathological images, an image in disease class could contain a large portion of healthy nuclei. This small between-class variance issue is another aspect we should consider.

Sparse representations have emerged as a powerful tool for a range of signal processing applications. Applications include compressed sensing [15], signal denoising, sparse signal recovery [16], image inpainting [17], image segmentation [18], and more recently, signal classification. In such representations, most of signals can be expressed by a linear combination of few bases taken from a “dictionary”. Based on this theory, a sparse representation-based classifier (SRC) [1] was initially developed for robust face recognition, and thereafter adapted to numerous signal/image classification problems, ranging from medical image classification [19, 2, 20], hyperspectral image classification [21, 22, 23], synthetic aperture radar (SAR) image classification [24], recaptured image recognition [25], video anomaly detection [26], and several others [27, 28, 29, 30, 31, 32].

The success of SRC-related methods firstly comes from the fact that sparsity models are often robust to noise and occlusions [1]. In addition, the insufficient training problem can be mitigated by incorporating prior knowledge of signals into the optimization problem during the inference process as regularization terms or sparsity constraints that capture signal relationships [27, 33, 20]. It has been shown that learning a dictionary from the training samples instead of using all of them as a dictionary can further enhance the performance of SRC.

This dissertation proposes novel structured sparsity frameworks for different classification applications. These applications include but are not limited to disease diagnosis and cancer detection, ultra-wide band synthetic aperture radar signal classification, face identification, flower classification, and general object classification. The works in this dissertation have culminated into one software and two toolboxes which are widely used by peer researchers.

The next section of this chapter provides overview of sparse representation-based classification and fundamental dictionary learning methods. The last section will present contributions and organization of this dissertation.

Sparse Representation-based Classification

Sparse Representation-based Classification

A significant contribution to the development of algorithms for image classification that addressed some of aforementioned challenges up to some extent is a recent sparse representation-based classification (SRC) framework [1], which exploits the discriminative capability of sparse representation. Given a sufficiently diverse collection of training images from each class, any image from a specific class can be approximately represented as a linear combination of training images from the same class. Therefore, if we have training images of all classes and form a basis or dictionary based on that, any new and unseen test image has a sparse representation with respect to such over-complete dictionary. It is worth to mention that sparsity assumption holds due to the class-specific design of dictionaries as well as the assumption of the linear representation model.

Concretely, given a collection $\mathbf{D}\in\mathbb{R}^{d\times k_{i}}$ of objects from one class, in which each column of $\mathbf{D}$ is the vectorized version of one signal (image in this case), a new signal $\mathbf{y}\in\mathbb{R}^{d}$ from the same class can be approximately expressed as $\mathbf{y}\approx\mathbf{D}\mathbf{x}$. In addition, this assumption is not true when $\mathbf{D}$ comprising images coming from other class. Now suppose that we have $C$ classes of subject, and let $\mathbf{D}=[\mathbf{D}_{1},\mathbf{D}_{2},\dots,\mathbf{D}_{C}]$ be the set of original training samples, where $\mathbf{D}_{i}$ is the sub-set of training samples from class $i$. Denote by $\mathbf{y}$ a testing sample. The procedures of SRC are as follows:

Sparsely code $\mathbf{y}$ on $\mathbf{D}$ via $l_{1}$-norm minimization:

[TABLE]

where $\lambda$ is a scalar constant. 2. 2.

Do classification via:

[TABLE]

where $e_{i}=\|\mathbf{y}-\mathbf{D}_{i}\delta_{i}{(\hat{\mathbf{x}})}\|$ and $\delta_{i}(\hat{\mathbf{x}})$ is the part of $\hat{\mathbf{x}}$ associated with class $i$.

Although SRC scheme shows interesting results, the dictionary used in it may not be effective enough to represent the query images due to the uncertain and noisy information in the original training images. The coding complexity increases as more training data is involved in building the dictionary. In addition, using the original training samples as the dictionary could not fully exploit the discriminative information hidden in the training samples. On the other hand, using analytically designed off-the-shelf bases as dictionary (e.g. [34] uses Haar wavelets and Gabor wavelets as the dictionary) might be universal to all types of images but will not be effective enough for specific type of images such as face, digit and texture images. In fact, all the above mentioned problems of predefined dictionary can be addressed, at least to some extent, by learning properly a dictionary from the original training samples. The next subsections will review the Online Dictionary Learning Method for compact purpose and two well-known Discriminative Dictionary Learning methods used for classification purpose.

Online Dictionary Learning

When number of training images increases, concatenating all of them into one “fat” matrix $\mathbf{D}$ and then solving the sparse coding problem (85) would be a time-consuming task. Additionally, it would be extremely redundant if some images in one class look very similar. To address this problem, several dictionary learning methods have been proposed for the problem of reconstruction. Concretely, from the big training set $\mathbf{Y}$ of one class, the dictionary learning algorithm tries to find a comprehensive set of bases which have ability to sparsely represent all element in the set. This task could be done by solving the following optimization problem:

[TABLE]

where columns of $\mathbf{D}$ are constrained by $\left\|\mathbf{d}_{j}\right\|_{2}^{2}\leq 1$ to avoid trivial solutions. Problem (3) could be rewritten in the matrix form:

[TABLE]

where $\left\|\mathbf{X}\right\|_{1}$ is the sum of absolute values of all elements in $\mathbf{X}$.

Problem (4) is not simultaneously convex with respect to both $\mathbf{D}$ and $\mathbf{X}$ but convex with each variable if the other variable is fixed. The typical approach is as follows. First, suppose that $\mathbf{D}$ is fixed, then $\mathbf{X}$ could be found using LASSO[35]:

[TABLE]

Second, fixing $\mathbf{X}=\mathbf{X}^{*}$ and compute $\mathbf{D}$ by:

[TABLE]

Mairal et al.[36] propose an algorithm for computing $\mathbf{D}$ by updating column by column until convergence:

[TABLE]

where $\mathbf{F}(j,j)$ is the value of $\mathbf{F}$ at coordinate $(j,j)$ and $\mathbf{f}_{j}$ denotes the $j$-th column of $\mathbf{F}$.

The dictionary learned from this method has been shown to be comprehensive in terms of sparsely expressing in-class samples. Because ODL is trained using in-class samples only, it might or might not well present complementary samples. However, from classification view point, this fact could be useful in terms of discriminability. Several dictionary learning methods have been proposed for the classification purpose. Two well-known dictionary learning methods demonstrating impressive results in object recognition are LC-KSVD[37] and FDDL[38], which are also reviewed in the following sections of this chapter.

Label Consistent K-SVD (LC-KSVD)

Z. Jiang et al.[37] propose another dictionary learning method which learns a single over-complete dictionary and an optimal linear classifier simultaneously. It yields dictionaries so that feature points with the same class labels have similar sparse codes. Each dictionary item is chosen so that it can be associated with a particular label. It is claimed that the performance of the linear classifier depends on the discriminability of the input sparse codes $\mathbf{x}$. For obtaining discriminative sparse codes $\mathbf{x}$ with the learned $\mathbf{D}$, an objective function for dictionary construction is defined as:

[TABLE]

where $\alpha$ controls the relative contribution between reconstruction and label consistent regularization, and $\mathbf{Q}=[\mathbf{q}_{1},\mathbf{q}_{2},\dots,\mathbf{q}_{N}]\in\mathbb{R}^{d\times N}$ are the ’discriminative’ sparse codes of input signals $\mathbf{Y}$ for classification. They say that $\mathbf{q}_{i}=[\mathbf{q}_{i}^{1},\mathbf{q}_{i}^{2},\dots,\mathbf{q}_{i}^{N}]^{T}=[0,\dots,1,1,\dots,0]^{T}\in\mathbb{R}^{d}$ is a ’discriminative’ sparse code corresponding to an input signal $\mathbf{y}_{i}$, if the non-zero values of $\mathbf{q}_{i}$ occur at those indicates where the input signal $\mathbf{y}_{i}$ and the dictionary item $\mathbf{d}_{k}$ share the same label.

Classification scheme: After the desired dictionary $\hat{\mathbf{D}}$ has been learned, for a test image $\mathbf{y}$, they first compute its sparse representation $\mathbf{x}$ by solving the optimization problem:

[TABLE]

Then label of the image $\mathbf{y}$ is estimated by using the linear predictive classifier $\hat{\mathbf{W}}$:

[TABLE]

where $l\in\mathbb{R}^{m}$ is the class label vector.

Fisher Discrimination Dictionary Learning (FDDL)

M. Yang et al.[38] proposed another dictionary learning method to improve the pattern classification performance compared to SRC[1]. A structured dictionary, whose bases have correspondence to the class labels, is learned so that the reconstruction error after sparse coding can be used for pattern classification. Meanwhile, the Fisher discrimination criterion is imposed on the coding coefficients so that they have small within-class scatter but big between-class scatter.

Specifically, the “total” dictionary $\mathbf{D}$ consisting of $c$ class-specific dictionaries $\mathbf{D}=[\mathbf{D}_{1},\mathbf{D}_{2},\dots,\mathbf{D}_{c}]$ is learned via the following optimization problem:

[TABLE]

where:

• •

$\displaystyle r(\mathbf{Y}_{i},\mathbf{D},\mathbf{X}_{i})=\|\mathbf{Y}_{i}-\mathbf{D}\mathbf{X}_{i}\|_{F}^{2}+\|\mathbf{Y}_{i}-\mathbf{D}_{i}\mathbf{X}_{i}^{i}\|_{F}^{2}+\sum_{j=1,j\neq i}^{c}\|\mathbf{D}_{j}\mathbf{X}_{i}^{j}\|_{F}^{2}$: the discriminative fidelity term.

• •

$\mathbf{X}_{i}^{j}$: the coding coefficient of $\mathbf{Y}_{i}$ over the sub-dictionary $\mathbf{D}_{j}$.

• •

$f(\mathbf{X})=\textrm{trace}(S_{W}(\mathbf{X})-S_{B}(\mathbf{X}))+\eta\|\mathbf{X}\|_{F}^{2}$: the discriminative coefficient term.

More details about the optimization of FDDL and its classification scheme could be found at [38].

Dissertation Contributions and Organization

A snapshot of the main contributions of this dissertation is presented next. Publications related to the contribution in each chapter are also listed where applicable.

In Chapter 2, the primary contribution is the proposal of a new Discriminative Feature-oriented Dictionary Learning (DFDL) method for automatic feature discovery in histopathological images. This proposal mitigates the generally difficulty of feature extraction in histopathological images. Our discriminative framework learns dictionaries that emphasize inter-class differences while keeping intra-class differences small, resulting in enhanced classification performance. The design is based on solving a sparsity constrained optimization problem, for which we develop a tractable algorithmic solution.

Experimental validation of DFDL is carried out on three diverse histopathological datasets to show its broad applicability. The first dataset is courtesy of the Clarian Pathology Lab and Computer and Information Science Dept., Indiana University-Purdue University Indianapolis (IUPUI). The images acquired by the process described in [39] correspond to human Intraductal Breast Lesions (IBL). Two well-defined categories will be classified: Usual Ductal Hyperplasia (UDH)–benign, and Ductal Carcinoma In Situ (DCIS)–actionable. The second dataset contains images of brain cancer (glioblastoma or GBM) obtaind from The Cancer Genome Atlas (TCGA) [40] provided by the National Institute of Health, and will henceforth be referred as the TCGA dataset. For this dataset, we address the problem of detecting MicroVascular Proliferation (MVP) regions, which is an important indicator of a high grade glioma (HGG) [41]. The third dataset is provided by the Animal Diagnostics Lab (ADL), The Pennsylvania State University. It contains tissue images from three mammalian organs - kidney, lung and spleen. For each organ, images will be assigned into one of two categories–healthy or inflammatory. The samples of these three datasets are given in Figs. 2,  3, and 4, respectively. Extensive experimental results show that our method outperforms many competing methods, particularly in low training scenarios. In addition, Receiver Operating Characteristic (ROC) curves are provided that facilitate a trade-off between false alarm and miss rates. This work was done in collaboration with Prof. Ganesh Rao at the Department of Neurosurgery, and Prof. UK Arvind Rao at the Department of Bionformatics and Computational Biology, both at the University of Texas MD Anderson Cancer Center, Houston, TX, USA.

We derive the computational complexity of DFDL as well as competing dictionary learning methods in terms of approximate number of operations needed. We also report experimental running time of DFDL and three other dictionary learning methods. All results in the manuscript are reproducible via a user-friendly software111The software can be downloaded at http://signal.ee.psu.edu/dfdl.html. The software (MATLAB toolbox) is also provided with the hope of usage in future research and comparisons via peer researchers.

This material was presented at the 2015 IEEE International Symposium on Biomedical Imaging and appeared in the IEEE Transactions of Medical Imaging in March 2016.

Chapter 3 proposes a new low-rank shared dictionary learning framework (LRSDL) for automatically extracting both discriminative and shared bases in several widely used image datasets to enhance the classification performance of dictionary learning methods. Our framework simultaneously learns each class-dictionary per class to extract discriminative features and the shared features that all classes contain. For the shared part, we impose two intuitive constraints. First, the shared dictionary must have a low-rank structure. Otherwise, the shared dictionary may also expand to contain discriminative features. Second, we contend that the sparse coefficients corresponding to the shared dictionary should be almost similar. In other words, the contribution of the shared dictionary to reconstruct every signal should be close together. We will experimentally show that both of these constraints are crucial for the shared dictionary.

Significantly, new accurate and efficient algorithms for selected existing and proposed dictionary learning methods are proposed. We present three effective algorithms for dictionary learning: i) sparse coefficient update in FDDL [4] by using FISTA [42]. We address the main challenge in this algorithm – how to calculate the gradient of a complicated function effectively – by introducing a new simple function $\mathcal{M}(\bullet)$ on block matrices and a lemma to support the result. ii) Dictionary update in FDDL [4] by a simple ODL [36] procedure using $\mathcal{M}(\bullet)$ and another lemma. Because it is an extension of FDDL, the proposed LRSDL also benefits from the aforementioned efficient procedures. iii) Dictionary update in DLSI [43] by a simple ADMM [44] procedure which requires only one matrix inversion instead of several matrix inversions as originally proposed in [43]. We subsequently show the proposed algorithms have both performance and computational benefits.

We derive the computational complexity of numerous dictionary learning methods in terms of approximate number of operations (multiplications) needed. We also report complexities and experimental running time of aforementioned efficient algorithms and their original counterparts. Numerous sparse coding and dictionary learning algorithms in the manuscript are reproducible via a user-friendly toolbox. The toolbox includes implementations of SRC [1], ODL [36], LC-KSVD [37]222Source code for LC-KSVD is directly taken from the paper at:

http://www.umiacs.umd.edu/$\sim$zhuolin/projectlcksvd.html., efficient DLSI [43], efficient COPAR [3], efficient FDDL [4], $D^{2}L^{2}R^{2}$ [45] and the proposed LRSDL. The toolbox (a MATLAB version and a Python version) is provided333The toolbox can be downloaded at: https://github.com/tiepvupsu/DICTOL. with the hope of usage in future research and comparisons via peer researchers.

The material in Chapter 3 was presented at the 2016 IEEE International Conference on Image Processing and appeared in the IEEE Transactions on Image Processing in November 2017.

In Chapter 4, a framework for simultaneously denoising and classifying 2-D UWB SAR imagery is introduced. Subtle features from targets of interest are directly learned from their SAR imagery. The classification also exploits polarization diversity and consecutive aspect angle dependence information of targets. A generalized tensor discriminative dictionary learning (TensorDL) is also proposed when more training data involved. These dictionary learning frameworks are shown to be robust even with high levels of noise.

Additionally, a relative SRC framework (ShiftSRC) is proposed to deal with multi-look data. Low-frequency UWB SAR signals are often captured at different views of objects, depending on the movement of the radar carriers. These signals contain uniquely important information of consecutive views. With ShiftSRC, this information will be comprehensively exploited. Importantly, a solution to the ShiftSRC framework can be obtained by an elegant modification on the training dictionary, resulting in a tensor sparse coding problem, which is similar to a problem proposed in the last paragraph. This work was done under the mentorship of Dr. Lam Nguyen at the U.S. Army Research Laboratory, Adelphi, MD. All tensor sparsity algorithms in this chapter are reproducible via a user-friendly toolbox. The toolbox written in Matlab is provided444 https://github.com/tiepvupsu/tensorsparsity with the hope of usage in future research via peer researchers.

The material in this chapter was presented at the 2017 IEEE Radar Conference and is under review at the IEEE Transactions on Aerospace and Electronic Systems. On a related note, we employed a deep learning framework for simultaneously denoising and classifying UWB SAR imagery. This work is however not included in this dissertation. This work was recently presented as an invited talk at the 2018 IEEE Radar Conference.

In Chapter 5, the main contributions of this dissertation are summarized.

As a side note, apart from building discriminative models for signal classification, our research also solves other long-standing open problems in sparse signal and image processing. Using sparsity as a prior is tremendously interesting in a wide variety of applications; however, existing solutions to address this issue are sub-optimal and often fail to capture the intrinsic sparse structure of physical phenomenon. We have been trying to address a very fundamental question in this area of how to efficiently and effectively capture sparsity in natural signals by using the Spike and Slab priors. These priors have been of much recent interest in signal processing as a means of inducing sparsity in Bayesian inference. It is well-known that solving for the sparse coefficient vector to maximize these priors results in a hard non-convex and mixed integer programming problem. Most existing solutions to this optimization problem either involve simplifying assumptions/relaxations or are computationally expensive. We propose a new greedy and adaptive matching pursuit (AMP) algorithm to directly solve this hard problem. Essentially, in each step of the algorithm, the set of active elements would be updated by either adding or removing one index, whichever results in better improvement. In addition, the intermediate steps of the algorithm are calculated via an inexpensive Cholesky decomposition which makes the algorithm much faster. Results on simulated data sets as well as real-world image recovery challenges confirm the benefits of the proposed AMP, particularly in providing a superior cost-quality trade-off over existing alternatives. These findings recently appeared in 2017 at the IEEE International Conference on Acoustics, Speech, and Signal Processing (ICASSP), and was nominated as a Finalist for the Best Student Paper Award. This work is nevertheless not included in this dissertation.

Discriminative Feature-oriented Dictionary Learning for Histopathological Image Classification

Introduction

Automated histopathological image analysis has recently become a significant research problem in medical imaging and there is an increasing need for developing quantitative image analysis methods as a complement to the effort of pathologists in diagnosis process. Consequently, an emerging class of problems in medical imaging focuses on the the development of computerized frameworks to classify histopathological images [46, 47, 20, 48, 41]. These advanced image analysis methods have been developed with three main purposes of (i) relieving the workload on pathologists by sieving out obviously diseased and also healthy cases, which allows specialists to spend more time on more sophisticated cases; (ii) reducing inter-expert variability; and (iii) understanding the underlying reasons for a specific diagnosis that pathologists might not realize.

In the diagnosis process, pathologists often look for problem-specific visual cues, or features, in histopathological images in order to categorize a tissue image as one of the possible categories. These features might come from the distinguishable characteristics of cells or nuclei, for example, size, shape or texture[39, 46]. They could also come from spatially related structures of cells [41, 49, 50, 20]. In some cancer grading problems, features might include the presence of particular regions [41, 51]. Consequently, different customized feature extraction techniques for a variety of problems have been developed based on these observed features[52, 53, 54, 55, 56]. Morphological image features have been utilized in medical image segmentation[57] for detection of vessel-like patterns. Wavelet features and histograms are also a popular choice of features for medical imaging[58, 59]. Graph-based features such as Delaunay triangulation, Vonoroi diagram, minimum spanning tree[49], query graphs[60] have been also used to exploit spatial structures. Orlov et al.[52, 53] have proposed a multi-purpose framework that collects texture information, image statistics and transforms domain coefficients to be set of features. For classification purposes, these features are combined with powerful classifiers such as neural networks or support vector machines (SVMs). Gurcan et al.[46] provided detailed discussion of feature and classifier selection for histopathological analysis.

Sparse representation frameworks have also been proposed for medical applications recently [20, 48, 61]. Specifically, Srinivas et al.[47, 20] presented a multi-channel histopathological image as a sparse linear combination of training examples under channel-wise constraints and proposed a residual-based classification technique. Yu et al.[62] proposed a method for cervigram segmentation based on sparsity and group clustering priors. Song et al.[63, 64] proposed a locality-constrained and a large-margin representation method for medical image classification. In addition, Parvin et al.[48] combined a dictionary learning framework with an autoencoder to learn sparse features for classification. Chang et al.[65] extended this work by adding a spatial pyramid matching to enhance the performance.

Challenges and Motivation

While histopathological analysis shares some traits with other image classification problems, there are also principally distinct challenges specific to histopathology. The central challenge comes from the geometric richness of tissue images, resulting in the difficulty of obtaining reliable discriminative features for classification. Tissues from different organs have structural and morphological diversity which often leads to highly customized feature extraction solutions for each problem and hence the techniques lack broad applicability.

Being mindful of the aforementioned challenges, we design via optimization, a discriminative dictionary for each class by imposing sparsity constraints that minimizes intra-class differences, while simultaneously emphasizing inter-class differences. On one hand, small intra-class differences encourage the comprehensibility of the set of learned bases, which has ability of representing in-class samples with only few bases (intra class sparsity). This encouragement forces the model to find the representative bases in that class. On the other hand, large inter-class differences prevent bases of a class from sparsely representing samples from other classes. Concretely, given a dictionary from a particular class $\mathbf{D}$ with $k$ bases and a certain sparsity level $L\ll k$, we define an $L$-subspace of $\mathbf{D}$ as a span of a subset of $L$ bases from $\mathbf{D}$. Our proposed Discriminative Feature-oriented Dictionary Learning (DFDL) aims to build dictionaries with this key property: any sample from a class is reasonably close to an $L$-subspace of the associated dictionary while a complementary sample is far from any $L$-subspace of that dictionary. Illustration of the proposed idea is shown in Fig. 5.

Contributions

Notation

The vectorization of a small block (or patch)555In our work, a training vector is obtained by vectorizing all three RGB channels followed by concatenating them together to have a long vector. extracted from an image is denoted as a column vector $\mathbf{y}\in\mathbb{R}^{d}$ which will be referred as a sample. In a classification problem where we have $c$ different categories, collection of all data samples from class $i$ ($i$ can vary between $1$ to $c$) forms the matrix $\mathbf{X}_{i}\in\mathbb{R}^{d\times N_{i}}$ and let \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{i}\in\mathbb{R}^{d\times\bar{N}_{i}} be the matrix containing all complementary data samples i.e. those that are not in class $i$. We denote by $\mathbf{D}_{i}\in\mathbb{R}^{d\times k_{i}}$ the dictionary of class $i$ that is desired to be learned through our DFDL method.

For a vector $\mathbf{x}\in\mathbb{R}^{k}$, we denote by $\|\mathbf{x}\|_{0}$ the number of its non-zero elements. The sparsity constraint of $\mathbf{x}$ can be formulated as $\|\mathbf{x}\|_{0}\leq L$. For a matrix $\mathbf{X}$ , $\|\mathbf{X}\|_{0}\leq L$ means that each column of $\mathbf{X}$ has no more than $L$ non-zero elements.

Discriminative Feature-oriented Dictionary Learning

We aim to build class-specific dictionaries $\mathbf{D}_{i}$ such that each $\mathbf{D}_{i}$ can sparsely represent samples from class $i$ but is poorly capable of representing its complementary samples with small number of bases. Concretely, for the learned dictionaries we need:

[TABLE]

where $L_{i}$ controls the sparsity level. These two sets of conditions could be simplified in the matrix form:

[TABLE]

The averaging operations $\left(\displaystyle\frac{1}{N_{i}}\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }\displaystyle\frac{1}{\bar{N_{i}}}\right)$ are taken here for avoiding the case where the largeness of inter-class differences is solely resulting from $\bar{N_{i}}\gg N_{i}$.

For simplicity, from now on, we consider only one class and drop the class index in each notion, i.e., using \mathbf{Y},\mathbf{D},\mathbf{X},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}},N,\bar{N},L instead of \mathbf{Y}_{i},\mathbf{D}_{i},\mathbf{X}_{i},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{i},N_{i},\bar{N}_{i} and $L_{i}$. Based on the argument above, we formulate the optimization problem for each dictionary:

[TABLE]

where $\rho$ is a positive regularization parameter. The first term in the above optimization problem encourages intra-class differences to be small, while the second term, with minus sign, emphasizes inter-class differences. By solving the above problem, we can jointly find the appropriate dictionaries as we desire in (15) and (16).

How to choose $L$: The sparsity level $L$ for classes might be different. For one class, if $L$ is too small, the dictionary might not appropriately express in-class samples, while if it is too large, the dictionary might be able to represent complementary samples as well. In both cases, the classifier might fail to determine identity of one new test sample. We propose a method for estimating $L$ as follows. First, a dictionary is learned using ODL[36] using in-class samples $\mathbf{Y}$ only:

[TABLE]

where $\lambda$ is a positive regularization parameter controlling the sparsity level. Note that the same $\lambda$ can still lead to different $L$ for different classes, depending on the intra-class variablity of each class. Without prior knowledge of those variablities, we choose the same $\lambda$ for every class. After $\mathbf{D}^{0}$ and $\mathbf{X}^{0}$ have been computed, $\mathbf{D}^{0}$ could be utilized as a warm initialization of $\mathbf{D}$ in our algorithm, $\mathbf{X}^{0}$ could be used to estimate the sparsity level $L$:

[TABLE]

Classification scheme: In the same manner with SRC [1], a new patch $\mathbf{y}$ is classified as follows. Firstly, the sparse codes $\hat{\mathbf{x}}$ are calculated via $l_{1}$-norm minimization:

[TABLE]

where $\mathbf{D}_{total}=[\mathbf{D}_{1},\mathbf{D}_{2},\dots,\mathbf{D}_{c}]$ is the collection of all dictionaries and $\gamma$ is a scalar constant. Secondly, the identity of $\mathbf{y}$ is determined as: $\displaystyle\arg\min_{i\in\{1,\dots,c\}}\{r_{i}(\mathbf{y})\}$ where

[TABLE]

and $\delta_{i}(\hat{\mathbf{x}})$ is part of $\hat{\mathbf{x}}$ associated with class $i$.

Proposed solution

We use an iterative method to find the optimal solution for the problem in (17). Specifically, the process is iterative by fixing $\mathbf{D}$ while optimizing $\mathbf{X}$ and  \overline{\hbox{}}$$\mathbf{X} and vice versa.

In the sparse coding step, with fixed $\mathbf{D}$, optimal sparse codes \mathbf{X}^{*},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}^{*} can be found by solving:

[TABLE]

With the same dictionary $\mathbf{D}$, these two sparse coding problems can be combined into the following one:

[TABLE]

with \hat{\mathbf{Y}}=[\mathbf{Y},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{Y}}}] being the matrix of all training samples and \hat{\mathbf{X}}=[\mathbf{X},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}]. This sparse coding problem can be solved effectively by OMP[66] using SPAMS toolbox[67].

For the bases update stage, $\mathbf{D}^{*}$ is found by solving:

[TABLE]

We have used the equation $\|\mathbf{M}\|_{F}^{2}=\textrm{trace}(\mathbf{M}\mathbf{M}^{T})$ for any matrix $\mathbf{M}$ to derive (24) from (23) and denoted:

[TABLE]

The objective function in (24) is very similar to the objective function in the dictionary update stage problem in [36] except that it is not guaranteed to be convex. It is convex if and only if $\mathbf{F}$ is positive semidefinite. For the discriminative dictionary learning problem, the symmetric matrix $\mathbf{F}$ is not guaranteed to be positive semidefinite, even all of its eigenvalues are real. In the worst case, where $\mathbf{F}$ is negative semidefinite, the objective function in (24) becomes concave; if we apply the same dictionary update algorithm as in [36], we will obtain its maximum solution instead of the minimum.

To deal with this situation, we propose a technique which convexifies the objective function based on the following observation.

If we look back to the main optimization problem stated in (17):

[TABLE]

we can see that if $\mathbf{D}=\left[\begin{matrix}\mathbf{d}_{1}&\mathbf{d}_{2}&\dots&\mathbf{d}_{k}\end{matrix}\right]$ is an optimal solution, then $\displaystyle\mathbf{D}=\left[\begin{matrix}\frac{\mathbf{d}_{1}}{a_{1}}&\frac{\mathbf{d}_{2}}{a_{2}}&\dots&\frac{\mathbf{d}_{k}}{a_{k}}\end{matrix}\right]$ is also an optimal solution as we multiply $j$-th rows of optimal $\mathbf{X}$ and  \overline{\hbox{}}$$\mathbf{X} by $a_{j}$, where $a_{j},j=1,2,\dots,k,$ are arbitrary nonzero scalars. Consequently, we can introduce constraints: $\|\mathbf{d}_{i}\|_{2}^{2}=1,j=1,2,\dots,k$, without affecting optimal value of (24). With these constraints, $\textrm{trace}(\mathbf{D}\lambda_{\min}(\mathbf{F})\mathbf{I}_{k}\mathbf{D}^{T})=\lambda_{\min}(\mathbf{F})\textrm{trace}(\mathbf{D}^{T}\mathbf{D})=\lambda_{\min}(\mathbf{F})\sum_{i=1}^{k}\mathbf{d}_{i}^{T}\mathbf{d}_{i}=k\lambda_{\min}(\mathbf{F})$, where $\lambda_{\min}(\mathbf{F})$ is the minimum eigenvalue of $\mathbf{F}$ and $\mathbf{I}_{k}$ denotes the identity matrix, is a constant. Substracting this constant from the objective function will not change the optimal solution to (24).

Essentially, the following problem in (26) is equivalent to (24):

[TABLE]

[TABLE]

The matrix $\hat{\mathbf{F}}=\mathbf{F}-\lambda_{\min}(\mathbf{F})\mathbf{I}_{k}$ is guaranteed to be positive semidefinite since all of its eignenvalues now are nonnegative, and hence the objective function in (26) is convex. Now, this optimization problem is very similar to the dictionary update problem in [36]. Then, $\mathbf{D}^{*}$ could be updated by the following iterations until convergence:

[TABLE]

where $\hat{\mathbf{F}}_{j,j}$ is the value of $\hat{\mathbf{F}}$ at coordinate $(j,j)$ and $\hat{\mathbf{f}}_{j}$ denotes the $j$-th column of $\hat{\mathbf{F}}$.

Our DFDL algorithm is summarized in Algorithm 1.

Overall classification procedures for three datasets

In this section, we propose a DFDL-based procedure for classifying images in three datasets.

IBL and ADL datasets

The key idea in this procedure is that a healthy tissue image largely consists of healthy patches which cover a dominant portion of the tissue. This procedure is shown in Fig. 6 and consists of the following three steps:

Step 1: Training DFDL bases for each class. From labeled training images, training patches are randomly extracted (they might be overlapping). The size of these patches is picked based on pathologist input and/or chosen by cross validation[68]. After we have a set of healthy patches and a set of diseased patches for training, class-specific DFDL dictionaries and the associated classifier are trained by using Algorithm 1.

Step 2: Learning a threshold $\theta$ for proportion of healthy patches in one healthy image. Labeled training images are now divided into non-overlapping patches. Each of these patches is then classified using the DFDL classifier as described in Eq. (20) and (21). The main purpose of this step is to find the threshold $\theta$ such that healthy images have proportion of healthy patches greater or equal to $\theta$ and diseased ones have proportion of diseased patches less than $\theta$. We can consider the proportion of healthy patches in one training image as its one-dimension feature. This feature is then put into a simple SVM to learn the threshold $\theta.$

Step 3: Classifying test images. For an unseen test image, we calculate the proportion $\tau$ of healthy patches in the same way described in Step 2. Now, the identity of the image is determined by comparing the proportion $\tau$ to $\theta$. It is categorized as healthy (diseased) if $\tau\geq(<)\theta$. The procedure readily generalizes to multi-class problems.

MVP detection problem in TCGA dataset

As described earlier, MicroVascular Proliferation (MVP) is the presence of blood vessels in a tissue and it is an important indicator of a high-grade tumor in brain glioma. Essentially presence of one such region in the tissue image indicates the high-grade tumor. Detection of such regions in TCGA dataset is an inherently hard problem and unlike classifying images in IBL and ADL datasets which are distinguishable by researching small regions, it requires more effort and investigation on larger connected regions. This is due to the fact that an MVP region may significantly vary in size and is usually surrounded by tumor cells which are actually benign or low grade. In addition, an MVP region is characterized by the presence of enlarged vessels in the tissue with different color shading and thick layers of cell rings inside the vessel (see Fig. 3). We define a patch as MVP if it lies entirely within an MVP region and as Not MVP otherwise. We also define a region as Not MVP if it does not contain any MVP patch. The procedure consists of two steps:

Step 1: Training phase. From training data, MVP regions and Not MVP regions are manually extracted. Note that while MVP regions come from MVP images only, Not MVP regions might appear in all images. From these extracted regions, DFDL dictionaries are obtained in the same way as in step 1 of IBL/ADL classification procedure described in section IBL and ADL datasets and Fig. 6.

**Step 2: ** MVP detection phase: A new unknown image is decomposed into non-overlapping patches. These patches are then classified using DFDL model learned before. After this step, we have a collection of patches classified as MVP. A region with large number of connected classified-as-MVP patches could be considered as an MVP region. If the final image does not contain any MVP region, we categorize the image as a Not MVP; otherwise, it is classified as MVP. The definition of connected regions contains a parameter $m$, which is the number of connected patches. Depending on $m$, positive patches might or might not appear in the final step. Specifically, if $m$ is small, false positives tend to be determined as MVP patches; if $m$ is large, true positives are highly likely eliminated. To determine $m$, we vary it from 1 to 20 and compute its ROC curve for training images and then simply pick the point which is closest to the origin and find the optimal $m$. This procedure is visualized in Fig. 7.

Validation and Experimental Results

In this section, we present the experimental results of applying DFDL to three diverse histopathological image datasets and compare our results with different competing methods:

$\bullet$ WND-CHARM[52, 53] in conjunction with SVM: this method combines state-of-the-art feature extraction and classification methods. We use the collection of features from WND-CHARM, which is known to be a powerful toolkit of features for medical images. While the original paper used weighted nearest neighbor as a classifier, we use a more powerful classifier (SVM [69]) to further enhance classification accuracy. We pick the most relevant features for histopathology[46], including but not limited to (color channel-wise) histogram information, image statistics, morphological features and wavelet coefficients from each color channel. The source code for WND-CHARM is made available by the National Institute of Health online at http://ome.grc.nia.nih.gov/.

$\bullet$ SRC[1]: We apply SRC on the vectorization of the luminance channel of the histopathological images, as proposed initially for face recognition and applied widely thereafter.

$\bullet$ SHIRC[20]: Srinivas et al.[47, 20] presented a simultaneous sparsity model for multi-channel histopathology image representation and classification which extends the standard SRC[1] approach by designing three color dictionaries corresponding to the RGB channels. The MATLAB code for the algorithms is posted online at: http://signal.ee.psu.edu/histimg.html.

$\bullet$ LC-KSVD[37] and FDDL[38]: These are two well-known dictionary learning methods which were applied to object recognition such as face, digit, gender, vehicle, animal, etc, but to our knowledge, have not been applied to histopathological image classification. To obtain a fair comparison, dictionaries are learned on the same training patches. Classification is then carried out using the learned dictionaries on non-overlapping patches in the same way described in Section Overall classification procedures for three datasets.

$\bullet$ Nayak’s: In recent relevant work, Nayak et al.[48] proposed a patch-based method to solve the problem of classification of tumor histopathology via sparse feature learning. The feature vectors are then fed into SVM to find the class label of each patch.

Experimental Set-Up: Image Datasets

IBL dataset: Each image contains a number of regions of interest (RoIs), and we have chosen a total of 120 images (RoIs), consisting of a randomly selected set of 20 images for training and the remaining 100 RoIs for test. Images are downsampled for computational purposes such that size of a cell is around 20-by-20 (pixels). Examples of images from this dataset are shown in Fig. 2. Experiments in section Validation of Central Idea: Visualization of Discovered Features below are conducted with 10 training images per class, 10000 patches of size 20-by-20 for training per class, $k=500$ bases for each dictionary, $\lambda=0.1$ and $\rho=0.001$. These parameters are chosen using cross-validation [68].

ADL dataset: This dataset contains bovine histopathology images from three sub-datasets of kidney, lung and spleen. Each sub-dataset consists of images of size $4000\times 3000$ pixels from two classes: healthy and inflammatory. Each class has around 150 images from which 40 images are chosen for training, the remaining ones are used for testing. Number of training patches, bases, $\lambda$ and $\rho$ are the same as in the IBL dataset. The classification procedure for IBL and ADL datasets is described in Section IBL and ADL datasets.

TCGA dataset: We use a total of 190 images (RoIs) (resolution $3000\times 3000$) from the TCGA, in which 57 images contain MVP regions and 133 ones have no traces of MVP. From each class, 20 images are randomly selected for training. The classification procedure for this dataset is described in Section MVP detection problem in TCGA dataset.

Each tissue specimen in these datasets is fixed on a scanning bed and digitized using a digitizer at 40$\times$ magnification.

Validation of Central Idea: Visualization of Discovered Features

This section provides experimental validation of the central hypothesis of this chapter: by imposing sparsity constraint on forcing intra-class differences to be small, while simultaneously emphasizing inter-class differences, the class-specific bases obtained are discriminative.

Example bases obtained by different dictionary learning methods are visualized in Fig. 8. By visualizing these bases, we emphasize that our DFDL is able to look for discriminative visual features from which pathologists could understand the reasons behind diseases. In the spleen dataset for example, it is really difficult to realize the differences between two classes by human eyes. However, by looking at DFDL learned bases, we can see that the distribution of cells in two classes are different such that a larger number of cells appears in a normal patch. These differences may provide pathologists one visual cue to classify these images without advanced tools. Moreover, for IBL dataset, UDH bases visualize elongated cells with sharp edges while DCIS bases present more rounded cells with blurry boundaries, which is consistent with their descriptions in [20] and [39]; for ADL-Lung, we observe that a healthy lung is characterized by large clear openings of the alveoli, while in the inflamed lung, the alveoli are filled with bluish-purple inflammatory cells. This distinction is very clear in the bases learned from DFDL where white regions appear more in normal bases than in inflammatory bases and no such information can be deduced from LC-KSVD or FDDL bases. In comparison, FDDL fails to discover discriminative visual features that are interpretable and LC-KSVD learns bases with the inter-class differences being less significant than DFDL bases. Furthermore, these LC-KSVD bases do not present key properties of each class, especially in lung dataset.

To understand more about the significance of discriminative bases for classification, let us first go back to SRC [1]. For simplicity, let us consider a problem with two classes with corresponding dictionaries $\mathbf{D}_{1}$ and $\mathbf{D}_{2}$. The identity of a new patch $\mathbf{y}$, which, for instance, comes from class 1, is determined by equations (20) and (21). In order to obtain good results, we expect most of active coefficients to be present in $\delta_{1}(\hat{\mathbf{x}})$. For $\delta_{2}(\hat{\mathbf{x}})$, its non-zeros, if they exists should have small magnitude. Now, suppose that one basis, $\mathbf{d}_{1}$, in $\mathbf{D}_{1}$ looks very similar to another basis, $\mathbf{d}_{2}$, in $\mathbf{D}_{2}$. When doing sparse coding, if one patch in class 1 uses $\mathbf{d}_{1}$ for reconstruction, it is highly likely that a similar patch $\mathbf{y}$ in the same class uses $\mathbf{d}_{2}$ for reconstruction instead. This misusage may lead to the case $\|\mathbf{y}-\mathbf{D}_{1}\delta_{1}(\hat{\mathbf{x}})\|>\|\mathbf{y}-\mathbf{D}_{2}\delta_{2}(\hat{\mathbf{x}})\|$, resulting in a misclassified patch. For this reason, the more discriminative bases are, the better the performance.

To formally verify this argument, we do one experiment on one normal and one inflammatory image from lung dataset in which the differences of DFDL bases and LCKSVD bases are most significant. From these images, patches are extracted, then their sparse codes are calculated using two dictionaries formed by DFDL bases and LC-KSVD bases. Fig. 9 demonstrates our results. Note that the plots in Figs. 9c) and d) are corresponding to DFDL while those in Figs. 9e) and f) are for LC-KSVD. Most of active coefficients in Fig. 9c) are gathered on the left of the red line, and their values are also greater than values on the right. This means that $\mathbf{D}_{1}$ contributes more to reconstructing the lung-normal image in Fig. 9a) than $\mathbf{D}_{2}$ does. Similarly, most of active coefficients in Fig. 9d) locate on the right of the vertical line. This agrees with what we expect since the image in Fig. 9a) belongs to class 1 and the one in Fig. 9b) belongs to class 2. On the contrary, for LC-KSVD, active coefficients in Fig. 9f) are more uniformly distributed on both sides of the red line, which adversely affects classification. In Fig. 9e), although active coefficients are strongly concentrated to the left of the red line, this effect is even more pronounced with DFDL, i.e. in Fig. 9c).

Overall Classification Accuracy

To verify the performance of our idea, for IBL and ADL datasets, we present overall classification accuracies in the form of bar graphs in Fig. 10. It is evident that DFDL outperforms other methods in both datasets. Specifically, in IBL and ADL Lung, the overall classification accuracies of DFDL are over 97.75$\%$, the next best rates come from WND-CHARM (92.85$\%$ in IBL) and FDDL (91.56$\%$ in ADL-Lung), respectively, and much higher than those reported in [20] and our own previous results in [19]. It is noteworthy to mention here that the overall classification accuracy of MIL[39] applied to IBL is under 90$\%$. In addition, for ADL-Kidney and ADL-Spleen, our DFDL also provides the best result with accuracy rates being nearly 90$\%$ and over 92$\%$, respectively.

For the TCGA dataset, overall accuracy of competing methods are shown in Fig. 11, which reveals that DFDL performance is the second best, bettered only by LC-KSVD and by less than 0.67$\%$ (i.e. one more misclassified image for DFDL).

Complexity analysis

In this section, we compare the computational complexity for the proposed DFDL and competing dictionary learning methods: LC-KSVD[37], FDDL[38], and Nayak’s[48]. The complexity for each dictionary learning method is estimated as the (approximate) number of operations required by each method in learning the dictionary (see Appendix Evolution of structured sparsity for details). From Table 1, it is clear that the proposed DFDL is the least expensive computationally. Note further, that the final column of Table 1 shows actual run times of each of the methods. The parameters were as follows: $c=2$ (classes), $k=500$ (bases per class), $N=10,000$ (training patches per class), data dimension $d=1200$ (3 channels $\times 20\times 20$), sparsity level $L=30$. The run time numbers in the final column of Table 1 are in fact consistent with numbers provided in Table 2, which are calculated by plugging the above parameters into the second column of Table 1.

Statistical Results: Confusion Matrices and ROC Curves

Next, we present a more elaborate interpretation of classification performance in the form of confusion matrices and ROC curves. Each row of a confusion matrix refers to the actual class identity of test images and each column indicates the classifier output. Table 3, 4 and 5 show the mean confusion matrices for all of three dataset. In continuation of trends from Fig. 10, in Table 4, DFDL offers the best disease detection accuracy in almost all datasets for each organ, while maintaining high classification accuracy for healthy images.

Typically in medical image classification problems, pathologists desire algorithms that reduce the probability of miss (diseased images are misclassified as healthy ones) while also ensuring that the false alarm rate remains low. However, there is a trade-off between these two quantities, conveniently described using receiver operating characteristic (ROC) curves. Fig. 12 and Fig. 11 (right) show the ROC curves for all three datasets. The lowest curve (closest to the origin) has the best overall performance and the optimal operating point minimizes the sum of the miss and false alarm probabilities. It is evident that ROC curves for DFDL perform best in comparison to those of other state-of-the-art methods.

Remark: Note for ROC comparisons, we compare the different flavors of dictionary learning methods (the proposed DFDL, LC-KSVD, FDDL and Nayak’s), this is because as Table 5 shows, they are the most competitive methods. Note for the IBL and ADL datasets, $\theta$, as defined in Fig. 6, is changed from 0 to 1 to acquire the curves; whereas for the TCGA dataset, number of connected classified-as-MVP patches, $m$, is changed from 1 to 20 to obtain the curves. It is worth re-emphasizing that DFDL achieves these results even as its complexity is lower than competing methods.

Performance vs. size of training set

Real-world histopathological classification tasks must often contend with lack of availability of large training sets. To understand training dependence of the various techniques, we present a comparison of overall classification accuracy as a function of the training set size for the different methods. We also present a comparison of classification rates as a function of the number of training patches for different dictionary learning methods666Since WND-CHARM is applied in the whole image level, there is no result for it in comparison of training patches.. In Fig. 13, overall classification accuracy is reported for IBL and ADL datasets corresponding to five scenarios. It is readily apparent that DFDL exhibits the most graceful decline as training is reduced.

Performance vs. number of training bases

We now compare the behavior of each dictionary learning method as the number of bases in each dictionary varies from 200 to 600 (with patch size being fixed at $20\times 20$ pixels). Results reported in Fig. 14 confirm that DFDL again outperforms other methods. In general, overall accuracies of DFDL on different datasets remain high when we reduce number of training bases. Interpreted another way, these results illustrate that DFDL is fairly robust to changes in parameters, which is a highly desirable trait in practice.

Discussion and Conclusion

In this chapter, we address the histopathological image classification problem from a feature discovery and dictionary learning standpoint. This is a very important and challenging problem and the main challenge comes from the geometrical richness of tissue images, resulting in the difficulty of obtaining reliable discriminative features for classification. Therefore, developing a framework capable of capturing this structural richness and being able to discriminate between different types is investigated and to this end, we propose the DFDL method which learns discriminative features for histopathology images. Our work aims to produce a more versatile histopathological image classification system through the design of discriminative, class-specific dictionaries which is hence capable of automatic feature discovery using example training image samples.

Our DFDL algorithm learns these dictionaries by leveraging the idea of sparse representation of in-class and out-of-class samples. This idea leads to an optimization problem which encourages intra-class similarities and emphasizes the inter-class differences. Ultimately, the optimization in (24) is done by solving the proposed equivalent optimization problem using a convexifying trick. Similar to other dictionary learning (machine learning approaches in general), DFDL also requires a set of regularization parameters. Our DFDL requires only one parameter, $\rho$, in its training process which is chosen by cross validation[68] – plugging different sets of parameters into the problem and selecting one which gives the best performance on the validation set. In the context of application of DFDL to real-world histopathological image slides, there are quite a few other settings should be carefully chosen, such as patch size, tiling method, number of connected components in the MVP detection etc. Of more importance is the patch size to be picked for each dataset which is mostly determined by consultation with the medical expert in the specific problem under investigation and the type of features that we should be looking for. For simplicity we employ regular tiling; however, using prior domain knowledge this may be improved. For instance in the context of MVP detection, informed selection of patch locations using existing disease detection and localization methods such as [41] can be used to further improve the detection of disease.

Experiments are carried out on three diverse histopathological datasets to show the broad applicability of the proposed DFDL method. It is illustrated our method is competitive with or outperforms state of the art alternatives, particularly in the regime of realistic or limited training set size. It is also shown that with minimal parameter tuning and algorithmic changes, DFDL method can be easily applied on different problems with different natures which makes it a good candidate for automated medical diagnosis instead of using customized and problem specific frameworks for every single diagnosis task. We also make a software toolbox available to help deploy DFDL widely as a diagnostic tool in existing histopathological image analysis systems. Particular problems such as grading and detecting specific regions in histopathology may be investigated using our proposed techniques.

Fast Low-rank Shared Dictionary Learning for Image Classification

Introduction

The central idea in SRC is to represent a test sample (e.g. a face) as a linear combination of samples from the available training set. Sparsity manifests because most of non-zeros correspond to bases whose memberships are the same as the test sample. Therefore, in the ideal case, each object is expected to lie in its own class subspace and all class subspaces are non-overlapping. Concretely, given $C$ classes and a dictionary $\mathbf{D}=[\mathbf{D}_{1},\dots,\mathbf{D}_{C}]$ with $\mathbf{D}_{c}$ comprising training samples from class $c,c=1,\dots,C$, a new sample $\mathbf{y}$ from class $c$ can be represented as $\mathbf{y}\approx\mathbf{D}_{c}\mathbf{x}^{c}$. Consequently, if we express $\mathbf{y}$ using the dictionary $\mathbf{D}:\mathbf{y}\approx\mathbf{D}\mathbf{x}=\mathbf{D}_{1}\mathbf{x}^{1}+\dots+\mathbf{D}_{c}\mathbf{x}^{c}+\dots+\mathbf{D}_{C}\mathbf{x}^{C}$, then most of active elements of $\mathbf{x}$ should be located in $\mathbf{x}^{c}$ and hence, the coefficient vector ${\mathbf{x}}$ is expected to be sparse. In matrix form, let $\mathbf{Y}=[\mathbf{Y}_{1},\dots,\mathbf{Y}_{c},\dots,\mathbf{Y}_{C}]$ be the set of all samples where $\mathbf{Y}_{c}$ comprises those in class $c$, the coefficient matrix $\mathbf{X}$ would be sparse. In the ideal case, $\mathbf{X}$ is block diagonal (see Figure 15).

Closely Related work and Motivation

The assumption made by most discriminative dictionary learning methods, i.e. non-overlapping subspaces, is unrealistic in practice. Often objects from different classes share some common features, e.g. background in scene classification. This problem has been partially addressed by recent efforts, namely DLSI [43], COPAR [3], JDL [33] and CSDL [70]. However, DLSI does not explicitly learn shared features since they are still hidden in the sub-dictionaries. COPAR, JDL and CSDL explicitly learn a shared dictionary $\mathbf{D}_{0}$ but suffer from the following drawbacks. First, we contend that the subspace spanned by columns of the shared dictionary must have low rank. Otherwise, class-specific features may also get represented by the shared dictionary. In the worst case, the shared dictionary span may include all classes, greatly diminishing the classification ability. Second, the coefficients (in each column of the sparse coefficient matrix) corresponding to the shared dictionary should be similar. This implies that features are shared between training samples from different classes via the “shared dictionary”. In this chapter, we develop a new low-rank shared dictionary learning framework (LRSDL) which satisfies the aforementioned properties. Our framework is basically a generalized version of the well-known FDDL [38, 4] with the additional capability of capturing shared features, resulting in better performance. We also show practical merits of enforcing these constraints are significant.

The typical strategy in optimizing general dictionary learning problems is to alternatively solve their subproblems where sparse coefficients $\mathbf{X}$ are found while fixing dictionary $\mathbf{D}$ or vice versa. In discriminative dictionary learning models, both $\mathbf{X}$ and $\mathbf{D}$ matrices furthermore comprise of several small class-specific blocks constrained by complicated structures, usually resulting in high computational complexity. Traditionally, $\mathbf{X}$, and $\mathbf{D}$ are solved block-by-block until convergence. Particularly, each block $\mathbf{X}_{c}$ (or $\mathbf{D}_{c}$ in dictionary update ) is solved by again fixing all other blocks $\mathbf{X}_{i},i\neq c$ (or $\mathbf{D}_{i},i\neq c$). Although this greedy process leads to a simple algorithm, it not only produces inaccurate solutions but also requires huge computation. In this chapter, we aim to mitigate these drawbacks by proposing efficient and accurate algorithms which allows to directly solve $\mathbf{X}$ and $\mathbf{D}$ in two fundamental discriminative dictionary learning methods: FDDL [4] and DLSI [43]. These algorithms can also be applied to speed-up our proposed LRSDL, COPAR [3], $D^{2}L^{2}R^{2}$ [45] and other related works.

Discriminative dictionary learning framework

Notation

In addition to notation stated in the Introduction, let $\mathbf{D}_{0}$ be the shared dictionary, $\mathbf{I}$ be the identity matrix with dimension inferred from context. For $c=1,\dots,C$; $i=0,1,\dots,C$, suppose that $\mathbf{Y}_{c}\in\mathbb{R}^{d\times n_{c}}$ and $\mathbf{Y}\in\mathbb{R}^{d\times N}$ with $N=\sum_{c=1}^{C}n_{c}$; $\mathbf{D}_{i}\in\mathbb{R}^{d\times k_{i}}$, $\mathbf{D}\in\mathbb{R}^{d\times K}$ with $K=\sum_{c=1}^{C}k_{c}$; and $\mathbf{X}\in\mathbb{R}^{K\times N}$. Denote by $\mathbf{X}^{i}$ the sparse coefficient of $\mathbf{Y}$ on $\mathbf{D}_{i}$, by $\mathbf{X}_{c}\in\mathbb{R}^{{K}\times N_{c}}$ the sparse coefficient of $\mathbf{Y}_{c}$ on ${\mathbf{D}}$, by $\mathbf{X}_{c}^{i}$ the sparse coefficient of $\mathbf{Y}_{c}$ on $\mathbf{D}_{i}$. Let \hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}}=\left[\begin{matrix}\mathbf{D}&\mathbf{D}_{0}\end{matrix}\right] be the total dictionary, \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}=[\mathbf{X}^{T},(\mathbf{X}^{0})^{T}]^{T} and \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{c}=[(\mathbf{X}_{c})^{T},(\mathbf{X}^{0}_{c})^{T}]^{T}. For every dictionary learning problem, we implicitly constrain each basis to have its Euclidean norm no greater than 1. These variables are visualized in Figure 16a).

Let $\mathbf{m},\mathbf{m}^{0},$ and $\mathbf{m}_{c}$ be the mean of $\mathbf{X},\mathbf{X}^{0},$ and $\mathbf{X}_{c}$ columns, respectively. Given a matrix $\mathbf{A}$ and a natural number $n$, define $\mu(\mathbf{A},n)$ as a matrix with $n$ same columns, each column being the mean vector of all columns of $\mathbf{A}$. If $n$ is ignored, we implicitly set $n$ as the number of columns of $\mathbf{A}$. Let $\mathbf{M}_{c}=\mu(\mathbf{X}_{c}),\mathbf{M}^{0}=\mu(\mathbf{X}^{0})$, and $\mathbf{M}=\mu(\mathbf{X},n)$ be the mean matrices. The number of columns $n$ depends on context, e.g. by writing $\mathbf{M}_{c}-\mathbf{M}$, we mean that $n=n_{c}$ . The ‘mean vectors’ are illustrated in Figure 16c).

Given a function $f(A,B)$ with $A$ and $B$ being two sets of variables, define $f_{A}(B)=f(A,B)$ as a function of $B$ when the set of variables $A$ is fixed. Greek letters ($\lambda,\lambda_{1},\lambda_{2},\eta$) represent positive regularization parameters. Given a block matrix $\mathbf{A}$, define a function $\mathcal{M}(\mathbf{A})$ as follows:

[TABLE]

That is, $\mathcal{M}(\mathbf{A})$ doubles diagonal blocks of $\mathbf{A}$. The row and column partitions of $\mathbf{A}$ are inferred from context. $\mathcal{M}(\mathbf{A})$ is a computationally inexpensive function of $\mathbf{A}$ and will be widely used in our LRSDL algorithm and the toolbox.

We also recall here the FISTA algorithm [42] for solving the family of problems:

[TABLE]

where $h(\mathbf{X})$ is convex, continuously differentiable with Lipschitz continuous gradient. FISTA is an iterative method which requires to calculate gradient of $h(\mathbf{X})$ at each iteration. In this chapter, we will focus on calculating the gradient of $h$.

Closely related work: Fisher discrimination dictionary learning (FDDL)

FDDL [38] has been used broadly as a technique for exploiting both structured dictionary and learning discriminative coefficient. Specifically, the discriminative dictionary $\mathbf{D}$ and the sparse coefficient matrix $\mathbf{X}$ are learned based on minimizing the following cost function:

[TABLE]

where $\displaystyle f_{\mathbf{Y}}(\mathbf{D},\mathbf{X})=\sum_{c=1}^{C}r_{\mathbf{Y}_{c}}(\mathbf{D},\mathbf{X}_{c})$ is the discriminative fidelity with:

$r_{\mathbf{Y}_{c}}(\mathbf{D},\mathbf{X}_{c})=\|\mathbf{Y}_{c}-\mathbf{D}\mathbf{X}_{c}\|_{F}^{2}+\|\mathbf{Y}_{c}-\mathbf{D}_{c}\mathbf{X}_{c}^{c}\|_{F}^{2}+\sum_{j\neq c}\|\mathbf{D}_{j}\mathbf{X}^{j}_{c}\|_{F}^{2},$

$g(\mathbf{X})=\sum_{c=1}^{C}(\|\mathbf{X}_{c}-\mathbf{M}_{c}\|_{F}^{2}-\|\mathbf{M}_{c}-\mathbf{M}\|_{F}^{2})+\|\mathbf{X}\|_{F}^{2}$ is the Fisher-based discriminative coefficient term, and the $l_{1}$-norm encouraging the sparsity of coefficients. The last term in $r_{\mathbf{Y}_{c}}(\mathbf{D},\mathbf{X}_{c})$ means that $\mathbf{D}_{j}$ has a small contribution to the representation of $\mathbf{Y}_{c}$ for all $j\neq c$. With the last term $\|\mathbf{X}\|_{F}^{2}$ in $g(\mathbf{X})$, the cost function becomes convex with respect to $\mathbf{X}$.

Proposed Low-rank shared dictionary learning (LRSDL)

The shared dictionary needs to satisfy the following properties:

1) Generativity: As the common part, the most important property of the shared dictionary is to represent samples from all classes [3, 33, 70]. In other words, it is expected that $\mathbf{Y}_{c}$ can be well represented by the collaboration of the particular dictionary $\mathbf{D}_{c}$ and the shared dictionary $\mathbf{D}_{0}$. Concretely, the discriminative fidelity term $f_{\mathbf{Y}}(\mathbf{D},\mathbf{X})$ in (31) can be extended to \hbox to0.0pt{\hskip 0.48611pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorf}}_{\mathbf{Y}}(\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}})=\sum_{c=1}^{C}\hskip 0.06319pt\hbox to0.0pt{\leavevmode\hbox{\set@colorr}\hss}{\hskip-0.06319pt\leavevmode\hbox{\set@color\overline{\hbox{}}}}_{\mathbf{Y}_{c}}(\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{c}) with \hskip 0.06319pt\hbox to0.0pt{\leavevmode\hbox{\set@colorr}\hss}{\hskip-0.06319pt\leavevmode\hbox{\set@color\overline{\hbox{}}}}_{\mathbf{Y}_{c}}(\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{c}) being defined as:

[TABLE]

Note that since \bar{r}_{\mathbf{Y}_{c}}(\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{c})=r_{\bar{\mathbf{Y}}_{c}}(\mathbf{D},\mathbf{X}_{c}) with \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{Y}}}_{c}=\mathbf{Y}_{c}-\mathbf{D}_{0}\mathbf{X}^{0}_{c} (see Figure 16b)), we have:

[TABLE]

with \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{Y}}}=\mathbf{Y}-\mathbf{D}_{0}\mathbf{X}^{0}.

This generativity property can also be seen in Figure 17a). In this figure, the intersection of different subspaces, each representing one class, is one subspace visualized by the light brown region. One class subspace, for instance class 1, can be well represented by the ideal shared atoms (dark brown triangles) and the corresponding class-specific atoms (red squares).

2) Low-rankness: The stated generativity property is only the necessary condition for a set of atoms to qualify a shared dictionary. Note that the set of atoms inside the shaded ellipse in Figure 17a) also satisfies the generativity property: along with the remaining red squares, these atoms well represent class 1 subspace; same can be observed for class 2. In the worst case, the set including all the atoms can also satisfy the generativity property, and in that undesirable case, there would be no discriminative features remaining in the class-specific dictionaries. Low-rankness is hence necessary to prevent the shared dictionary from absorbing discriminative atoms. The constraint is natural based on the observation that the subspace spanned by the shared dictionary has low dimension. Concretely, we use the nuclear norm regularization $\|\mathbf{D}_{0}\|_{*}$, which is the convex relaxation of $\text{rank}(\mathbf{D}_{0})$ [71], to force the shared dictionary to be low-rank. In contrast with our work, existing approaches that employ shared dictionaries, i.e. COPAR [3] and JDL [33], do not incorporate this crucial constraint.

3) Code similarity: In the classification step, a test sample $\mathbf{y}$ is decomposed into two parts: the part represented by the shared dictionary $\mathbf{D}_{0}\mathbf{x}^{0}$ and the part expressed by the remaining dictionary $\mathbf{D}\mathbf{x}$. Because $\mathbf{D}_{0}\mathbf{x}^{0}$ is not expected to contain class-specific features, it can be excluded before doing classification. The shared code $\mathbf{x}^{0}$ can be considered as the contribution of the shared dictionary to the representation of $\mathbf{y}$. Even if the shared dictionary already has low-rank, its contributions to each class might be different as illustrated in Figure 17b), the top row. In this case, the different contributions measured by $\mathbf{X}^{0}$ convey class-specific features, which we aim to avoid. Naturally, the regularization term $\|\mathbf{X}^{0}-\mathbf{M}^{0}\|$ is added to our proposed objective function to force each $\mathbf{x}^{0}$ to be close to the mean vector $\mathbf{m}^{0}$ of all $\mathbf{X}^{0}$.

With this constraint, the Fisher-based discriminative coefficient term $g(\mathbf{X})$ is extended to \hbox to0.0pt{\hskip 0.06422pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorg}}(\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}) defined as:

[TABLE]

Altogether, the cost function \hbox to0.0pt{\hskip 0.75346pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorJ}}_{\mathbf{Y}}(\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}) of our proposed LRSDL is:

[TABLE]

By minimizing this objective function, we can jointly find the class specific and shared dictionaries. Notice that if there is no shared dictionary $\mathbf{D}_{0}$ (by setting $k_{0}=0$), then \hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}} become $\mathbf{D},\mathbf{X}$, respectively, \hbox to0.0pt{\hskip 0.75346pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorJ}}_{\mathbf{Y}}(\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}) becomes $J_{\mathbf{Y}}(\mathbf{D},\mathbf{X})$ and LRSDL reduces to FDDL.

Classification scheme:

After the learning process, we obtain the total dictionary  \overline{\hbox{}}$$\mathbf{D} and mean vectors $\mathbf{m}_{c},\mathbf{m}^{0}$. For a new test sample $\mathbf{y}$, first we find its coefficient vector \hbox to0.0pt{\hskip 0.13889pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{x}}}=[\mathbf{x}^{T},(\mathbf{x}^{0})^{T}]^{T} with the sparsity constraint on  \overline{\hbox{}}$$\mathbf{x} and further encourage $\mathbf{x}^{0}$ to be close to $\mathbf{m}^{0}$:

[TABLE]

Using  \overline{\hbox{}}$$\mathbf{x} as calculated above, we extract the contribution of the shared dictionary to obtain \hbox to0.0pt{\hskip 0.13889pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{y}}}=\mathbf{y}-\mathbf{D}_{0}\mathbf{x}^{0}. The identity of $\mathbf{y}$ is determined by:

[TABLE]

where $w\in[0,1]$ is a preset weight for balancing the contribution of the two terms.

Efficient solutions for optimization problems

Before diving into minimizing the LRSDL objective function in (33), we first present efficient algorithms for minimizing the FDDL objective function in (31).

Efficient FDDL dictionary update

Recall that in [38], the dictionary update step is divided into $C$ subproblems, each updates one class-specific dictionary $\mathbf{D}_{c}$ while others fixed. This process is repeated until convergence. This approach is not only highly time consuming but also inaccurate. We will see this in a small example presented in Section Validation of efficient algorithms. We refer this original FDDL dictionary update as O-FDDL-D.

We propose here an efficient algorithm for updating dictionary called E-FDDL-D where the total dictionary $\mathbf{D}$ will be optimized when $\mathbf{X}$ is fixed, significantly reducing the computational cost.

Concretely, when we fix $\mathbf{X}$ in equation (31), the problem of solving $\mathbf{D}$ becomes:

[TABLE]

Therefore, $\mathbf{D}$ can be solved by using the following lemma.

Lemma 1**.**

The optimization problem (36) is equivalent to:

[TABLE]

where $\mathbf{E}=\mathbf{Y}\mathcal{M}(\mathbf{X}^{T})$ and $\mathbf{F}=\mathcal{M}(\mathbf{X}\mathbf{X}^{T})$.

Proof: See Appendix Proof of Lemma 1.

The problem (37) can be solved effectively by Online Dictionary Learning (ODL) method [36].

Efficient FDDL sparse coefficient update (E-FDDL-X)

When $\mathbf{D}$ is fixed, $\mathbf{X}$ will be found by solving:

[TABLE]

where $h(\mathbf{X})=\frac{1}{2}f_{\mathbf{Y},\mathbf{D}}(\mathbf{X})+\frac{\lambda_{2}}{2}g(\mathbf{X})$. The problem (38) has the form of equation (30), and can hence be solved by FISTA [42]. We need to calculate gradient of $f(\bullet)$ and $g(\bullet)$ with respect to $\mathbf{X}$.

Lemma 2**.**

Calculating gradient of $h(\mathbf{X})$ in equation (2)

[TABLE]

Then we obtain:

[TABLE]

Proof: See Appendix Proof of Lemma 2.

Since the proposed LRSDL is an extension of FDDL, we can also extend these two above algorithms to optimize LRSDL cost function as follows.

LRSDL dictionary update (LRSDL-D)

Returning to our proposed LRSDL problem, we need to find \hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}}=[\mathbf{D},\mathbf{D}_{0}] when  \overline{\hbox{}}$$\mathbf{X} is fixed. We propose a method to solve $\mathbf{D}$ and $\mathbf{D}_{0}$ separately.

For updating $\mathbf{D}$, recall the observation that \hbox to0.0pt{\hskip 0.48611pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorf}}_{\mathbf{Y}}(\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}})=f_{\hbox to0.0pt{\hskip 0.12498pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{Y}}}}(\mathbf{D},\mathbf{X}), with \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{Y}}}\triangleq\mathbf{Y}-\mathbf{D}_{0}\mathbf{X}^{0} (see equation (Proposed Low-rank shared dictionary learning (LRSDL))), and the E-FDDL-D presented in section Efficient FDDL dictionary update, we have:

[TABLE]

with \mathbf{E}=\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{Y}}}\mathcal{M}(\mathbf{X}^{T}) and $\mathbf{F}=\mathcal{M}(\mathbf{X}\mathbf{X}^{T}).$

For updating $\mathbf{D}_{0}$, we use the following lemma:

Lemma 3**.**

When $\mathbf{D},\mathbf{X}$ in (33) are fixed,

[TABLE]

where $\mathbf{V}=\mathbf{Y}-\frac{1}{2}\mathbf{D}\mathcal{M}(\mathbf{X})$.

Proof: See Appendix Proof of Lemma 3.

Based on the Lemma 3, $\mathbf{D}_{0}$ can be updated by solving:

[TABLE]

using the ADMM [44] method and the singular value thresholding algorithm [72]. The ADMM procedure is as follows. First, we choose a positive $\rho$, initialize $\mathbf{Z}=\mathbf{U}=\mathbf{D}_{0}$, then alternatively solve each of the following subproblems until convergence:

[TABLE]

where $\mathcal{D}$ is the shrinkage thresholding operator [72]. The optimization problem (45) can be solved by ODL [36]. Note that (46) and (48) are computationally inexpensive.

LRSDL sparse coefficients update (LRSDL-X)

In our preliminary work [73], we proposed a method for effectively solving $\mathbf{X}$ and $\mathbf{X}^{0}$ alternatively, now we combine both problems into one and find  \overline{\hbox{}}$$\mathbf{X} by solving the following optimization problem:

[TABLE]

where \displaystyle\hbox to0.0pt{\hskip 0.38078pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorh}}(\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}})=\frac{1}{2}\hbox to0.0pt{\hskip 0.48611pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorf}}_{\mathbf{Y},\overline{\mathbf{D}}}(\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}})+\frac{\lambda_{2}}{2}\hbox to0.0pt{\hskip 0.06422pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorg}}(\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}). We again solve this problem using FISTA [42] with the gradient of \displaystyle\hbox to0.0pt{\hskip 0.38078pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorh}}(\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}):

[TABLE]

For the upper term, by combining the observation

[TABLE]

and using equation, we obtain:

[TABLE]

For the lower term, by using Lemma 3, we have:

[TABLE]

[TABLE]

By combining these two terms, we can calculate (50).

Having \nabla\hbox to0.0pt{\hskip 0.38078pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorh}}(\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}) calculated, we can update  \overline{\hbox{}}$$\mathbf{X} by the FISTA algorithm [42] as given in Algorithm 1. Note that we need to compute a Lipschitz coefficient $L$ of \nabla\hbox to0.0pt{\hskip 0.38078pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@colorh}}(\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}). The overall algorithm of LRSDL is given in Algorithm 2.

Efficient solutions for other dictionary learning methods

We also propose here another efficient algorithm for updating dictionary in two other well-known dictionary learning methods: DLSI [43] and COPAR [3].

The cost function $J_{1}(\mathbf{D},\mathbf{X})$ in DLSI is defined as:

[TABLE]

Each class-specific dictionary $\mathbf{D}_{c}$ is updated by fixing others and solve:

[TABLE]

with $\mathbf{A}=\left[\begin{matrix}\mathbf{D}_{1},\dots,\mathbf{D}_{c-1},\mathbf{D}_{c+1},\dots,\mathbf{D}_{C}\end{matrix}\right]^{T}$.

The original solution for this problem, which will be referred as O-FDDL-D, updates each column $\mathbf{d}_{c,j}$ of $\mathbf{D}_{c}$ one by one based on the procedure:

[TABLE]

where $\mathbf{d}_{c,i}$ is the $i$-th column of $\mathbf{D}_{c}$ and $\mathbf{x}_{c}^{j}$ is the $j$-th row of $\mathbf{X}_{c}$. This algorithm is highly computational since it requires one matrix inversion for each of $k_{c}$ columns of $\mathbf{D}_{c}$. We propose one ADMM [44] procedure to update $\mathbf{D}_{c}$ which requires only one matrix inversion, which will be referred as E-DLSI-D. First, by letting $\mathbf{E}=\mathbf{Y}_{c}(\mathbf{X}^{c})^{T}$ and $\mathbf{F}=\mathbf{X}^{c}(\mathbf{X}^{c})^{T}$, we rewrite (56) in a more general form:

[TABLE]

In order to solve this problem, first, we choose a $\rho$, let $\mathbf{Z}=\mathbf{U}=\mathbf{D}_{c}$, then alternatively solve each of the following sub problems until convergence:

[TABLE]

This efficient algorithm requires only one matrix inversion. Later in this chapter, we will both theoretically and experimentally show that E-DLSI-D is much more efficient than O-DLSI-D [43]. Note that this algorithm can be beneficial for two subproblems of updating the common dictionary and the particular dictionary in COPAR [3] as well.

Complexity analysis

We compare the computational complexity for the efficient algorithms and their corresponding original algorithms. We also evaluate the total complexity of the proposed LRSDL and competing dictionary learning methods: DLSI [43], COPAR [3] and FDDL [38]. The complexity for each algorithm is estimated as the (approximate) number of multiplications required for one iteration (sparse code update and dictionary update). For simplicity, we assume: i) number of training samples, number of dictionary bases in each class (and the shared class) are the same, which means: $n_{c}=n,k_{i}=k$. ii) The number of bases in each dictionary is comparable to number of training samples per class and much less than the signal dimension, i.e. $k\approx n\ll d$. iii) Each iterative algorithm requires $q$ iterations to convergence. For consistency, we have changed notations in those methods by denoting $\mathbf{Y}$ as training sample and $\mathbf{X}$ as the sparse code.

In the following analysis, we use the fact that: i) if $\mathbf{A}\in\mathbb{R}^{m\times n},\mathbf{B}\in\mathbb{R}^{n\times p}$, then the matrix multiplication $\mathbf{A}\mathbf{B}$ has complexity $mnp$. ii) If $\mathbf{A}\in\mathbb{R}^{n\times n}$ is nonsingular, then the matrix inversion $\mathbf{A}^{-1}$ has complexity $n^{3}$. iii) The singular value decomposition of a matrix $\mathbf{A}\in\mathbb{R}^{p\times q}$, $p>q$, is assumed to have complexity $O(pq^{2})$.

Online Dictionary Learning (ODL)

We start with the well-known Online Dictionary Learning [36] whose cost function is:

[TABLE]

where $\mathbf{Y}\in\mathbb{R}^{d\times n},\mathbf{D}\in\mathbb{R}^{d\times k},\mathbf{X}\in\mathbb{R}^{k\times n}$. Most of dictionary learning methods find their solutions by alternatively solving one variable while fixing others. There are two subproblems:

Update $\mathbf{X}$ (ODL-X)

When the dictionary $\mathbf{D}$ is fixed, the sparse coefficient $\mathbf{X}$ is updated by solving the problem:

[TABLE]

using FISTA [42]. In each of $q$ iterations, the most computational task is to compute $\mathbf{D}^{T}\mathbf{D}\mathbf{X}-\mathbf{D}^{T}\mathbf{Y}$ where $\mathbf{D}^{T}\mathbf{D}$ and $\mathbf{D}^{T}\mathbf{Y}$ are precomputed with complexities $k^{2}d$ and $kdn$, respectively. The matrix multiplication $(\mathbf{D}^{T}\mathbf{D})\mathbf{X}$ has complexity $k^{2}n$. Then, the total complexity of ODL-X is:

[TABLE]

Update $\mathbf{D}$ (ODL-D)

After finding $\mathbf{X}$, the dictionary $\mathbf{D}$ will be updated by:

[TABLE]

subject to: $\|\mathbf{d}_{i}\|_{2}\leq 1$, with $\mathbf{E}=\mathbf{Y}\mathbf{X}^{T}$, and $\mathbf{F}=\mathbf{X}\mathbf{X}^{T}$.

Each column of $\mathbf{D}$ will be updated by fixing all others:

[TABLE]

where $\mathbf{d}_{i},\mathbf{e}_{i},\mathbf{f}_{i}$ are the $i-$th columns of $\mathbf{D},\mathbf{E},\mathbf{F}$ and $\mathbf{F}_{ii}$ is the $i-$th element in the diagonal of $\mathbf{F}$. The dominant computational task is to compute $\mathbf{D}\mathbf{f}_{i}$ which requires $dk$ operators. Since $\mathbf{D}$ has $k$ columns and the algorithm requires $q$ iterations, the complexity of ODL-D is $qdk^{2}$.

Dictionary learning with structured incoherence (DLSI)

DLSI [43] proposed a method to encourage the independence between bases of different classes by minimizing coherence between cross-class bases. The cost function $J_{1}(\mathbf{D},\mathbf{X})$ of DLSI is defined as (55).

Update $\mathbf{X}$ (DLSI-X)

In each iteration, the algorithm solves $C$ subproblems:

[TABLE]

with $\mathbf{Y}_{c}\in\mathbb{R}^{d\times n},\mathbf{D}_{c}\in\mathbb{R}^{d\times k}$, and $\mathbf{X}^{c}\in\mathbb{R}^{k\times n}$. Based on (66), the complexity of updating $\mathbf{X}$ ($C$ subproblems) is:

[TABLE]

Original update $\mathbf{D}$ (O-DLSI-D)

For updating $\mathbf{D}$, each sub-dictionary $\mathbf{D}_{c}$ is solved via (56). The main step in the algorithm is stated in (57) and (58). The dominant computational part is the matrix inversion which has complexity $d^{3}$. Matrix-vector multiplication and vector normalization can be ignored here. Since $\mathbf{D}_{c}$ has $k$ columns, and the algorithm requires $q$ iterations, the complexity of the O-DLSI-D algorithm is $Cqkd^{3}$.

Efficient update $\mathbf{D}$ (E-DLSI-D)

Main steps of the proposed algorithm are presented in equations (60)–(63) where (61) and (63) require much less computation compared to (60) and (62). The total (estimated) complexity of efficient $\mathbf{D}_{c}$ update is a summation of two terms: i) $q$ times ($q$ iterations) of ODL-D in (60). ii) One matrix inversion ($d^{3}$) and $q$ matrix multiplications in (62). Finally, the complexity of E-DLSI-D is:

[TABLE]

Total complexities of O-DLSI (the combination of DLSI-X and O-DLSI-D) and E-DLSI (the combination of DLSI-X and E-DLSI-D) are summarized in Table 7.

Separating the particularity and the commonality dictionary learning (COPAR)

Cost function

COPAR [3] is another dictionary learning method which also considers the shared dictionary (but without the low-rank constraint). By using the same notation as in LRSDL, we can rewrite the cost function of COPAR in the following form:

[TABLE]

where \displaystyle f_{1}(\mathbf{Y},\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}})=\sum_{c=1}^{C}r_{1}(\mathbf{Y}_{c},\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{c}) and r_{1}(\mathbf{Y}_{c},\hbox to0.0pt{\hskip 1.31944pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{D}}},\hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{c}) is defined as

[TABLE]

Update $\mathbf{X}$ (COPAR-X)

In sparse coefficient update step, COPAR [3] solve \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{X}}}_{c} one by one via one $l_{1}$-norm regularization problem:

[TABLE]

where $\tilde{\mathbf{Y}}\in\mathbb{R}^{\tilde{d}\times n},\tilde{\mathbf{D}}\in\mathbb{R}^{\tilde{d}\times\tilde{k}}$, $\tilde{\mathbf{X}}\in\mathbb{R}^{(\tilde{k}\times n}$, $\tilde{d}=2d+(C-1)k$ and $\tilde{k}=(C+1)k$ (details can be found in Section 3.1 of [3]). Following results in Section Update $\mathbf{X}$ (ODL-X) and supposing that $C\gg 1,q\gg 1,n\approx k\ll d$, the complexity of COPAR-X is:

[TABLE]

Update $\mathbf{D}$ (COPAR-D)

The COPAR dictionary update algorithm requires to solve $(C+1)$ problems of form (59). While O-COPAR-D uses the same method as O-DLSI-D (see equations (57-58)), the proposed E-COPAR-D takes advantages of E-DLSI-D presented in Section Efficient solutions for other dictionary learning methods. Therefore, the total complexity of O-COPAR-D is roughly $Cqkd^{3}$, while the total complexity of E-COPAR-D is roughly $C(q^{2}dk^{2}+d^{3}+qd^{2}k).$ Here we have supposed $C+1\approx C$ for large $C$.

Total complexities of O-COPAR (the combination of COPAR-X and O-COPAR-D) and E-COPAR (the combination of COPAR-X and E-COPAR-D) are summarized in Table 7.

Fisher discrimination dictionary learning (FDDL)

Original update $\mathbf{X}$ (O-FDDL-X)

Based on results reported in DFDL [2], the complexity of O-FDDL-X is roughly $C^{2}kn(d+qCk)+C^{3}dk^{2}=C^{2}k(dn+qCkn+Cdk)$.

Efficient update $\mathbf{X}$ (E-FDDL-X)

Based on section Efficient FDDL sparse coefficient update (E-FDDL-X), the complexity of E-FDDL-X mainly comes from equation (41). Recall that function $\mathcal{M}(\bullet)$ does not require much computation. The computation of $\mathbf{M}$ and ${\mathbf{M}_{c}}$ can also be neglected since each required calculation of one column, all other columns are the same. Then the total complexity of the algorithm E-FDDL-X is roughly:

[TABLE]

Original update $\mathbf{D}$ (O-FDDL-D)

The original dictionary update in FDDL is divided in to $C$ subproblems. In each subproblem, one dictionary $\mathbf{D}_{c}$ will be solved while all others are fixed via:

[TABLE]

where:

[TABLE]

When $d\gg k,C\gg 1$, complexity of updating $\mathbf{D}_{c}$ is:

[TABLE]

Then, complexity of O-FDDL-D is $Cdk(qk+C^{2}n)$.

Efficient update $\mathbf{D}$ (E-FDDL-D)

Based on Lemma 1, the complexity of E-FDDL-D is:

[TABLE]

Total complexities of O-FDDL and E-FDDL are summarized in Table 7.

LRSDL

Update $\mathbf{X},\mathbf{X}^{0}$

From (50), (52) and (54), in each iteration of updating  \overline{\hbox{}}$$\mathbf{X}, we need to compute:

[TABLE]

Therefore, the complexity of LRSDL-X is:

[TABLE]

which is similar to the complexity of E-FDDL-X. Recall that we have supposed number of classes $C\gg 1$.

Update $\mathbf{D}$

Compare to E-FDDL-D, LRSDL-D requires one more computation of \hbox to0.0pt{\hskip 1.25pt\leavevmode\hbox{\set@color\overline{\hbox{}}}\hss}{\leavevmode\hbox{\set@color\mathbf{Y}}}=\mathbf{Y}-\mathbf{D}_{0}\mathbf{X}^{0} (see section LRSDL dictionary update (LRSDL-D)). Then, the complexity of LRSDL-D is:

[TABLE]

which is similar to the complexity of E-FDDL-D.

Update $\mathbf{D}_{0}$

The algorithm of LRSDL-D0 is presented in section LRSDL dictionary update (LRSDL-D) with the main computation comes from (44), (45) and (47). The shrinkage thresholding operator in (47) requires one SVD and two matrix multiplications. The total complexity of LRSDL-D0 is:

[TABLE]

By combing (LABEL:eqn:complexity_lrsdl_X), (76) and (77), we obtain the total complexity of LRSDL, which is specified in the last row of Table 7.

Summary

Table 6 and Table 7 show final complexity analysis of each proposed efficient algorithm and their original counterparts. Table 7 compares LRSDL to other state-of-the-art methods. We pick a typical set of parameters with 100 classes, 20 training samples per class, 10 bases per sub-dictionary and shared dictionary, data dimension 500 and 50 iterations for each iterative method. Concretely, $C=100,\leavevmode\nobreak\ n=20,\leavevmode\nobreak\ k=10,\leavevmode\nobreak\ q=50,\leavevmode\nobreak\ d=500$. We also assume that in (77), $q_{2}=50$. Table 6 shows that all three proposed efficient algorithms require less computation than original versions with most significant improvements for speeding up DLSI-D. Table 7 demonstrates an interesting fact. LRSDL is the least expensive computationally when compared with other original dictionary learning algorithms, and only E-FDDL has lower complexity, which is to be expected since the FDDL cost function is a special case of the LRSDL cost function. COPAR is found to be the most expensive computationally.

Experimental results

Comparing methods and datasets

We present the experimental results of applying these methods to five diverse datasets: the Extended YaleB face dataset [74], the AR face dataset [75], the AR gender dataset, the Oxford Flower dataset [76], and two multi-class object category dataset – the Caltech 101 [77] and COIL-100 [78]. Example images from these datasets are shown in Figure 18. We compare our results with those using SRC [1] and other state-of-the-art dictionary learning methods: LC-KSVD [37], DLSI [43], FDDL [4], COPAR [3], $D^{2}L^{2}R^{2}$ [45], DLRD [79], JDL [33], and SRRS [80]. Regularization parameters in all methods are chosen using five-fold cross-validation [68]. For each experiment, we use 10 different randomly split training and test sets and report averaged results.

For two face datasets, feature descriptors are random faces, which are made by projecting face images onto a random vector using a random projection matrix. As in [81], the dimension of a random-face feature in the Extended YaleB is $d=504$, while the dimension in AR face is $d=540$. Samples of these two datasets are shown in Figure 18a) and b).

For the AR gender dataset, we first choose a non-occluded subset (14 images per person) from the AR face dataset, which consists of 50 males and 50 females, to conduct experiment of gender classification. Training images are taken from the first 25 males and 25 females, while test images comprises all samples from the remaining 25 males and 25 females. PCA was used to reduce the dimension of each image to 300. Samples of this dataset are shown in Figure 18c).

The Oxford Flower dataset is a collection of images of flowers drawn from 17 species with 80 images per class, totaling 1360 images. For feature extraction, based on the impressive results presented in [4], we choose the Frequent Local Histogram feature extractor [82] to obtain feature vectors of dimension 10,000. The test set consists of 20 images per class, the remaining 60 images per class are used for training. Samples of this dataset are shown in Figure 18d).

For the Caltech 101 dataset, we use a dense SIFT (DSIFT) descriptor. The DSIFT descriptor is extracted from $25\times 25$ patch which is densely sampled on a dense grid with 8 pixels. We then extract the sparse coding spatial pyramid matching (ScSPM) feature [83], which is the concatenation of vectors pooled from words of the extracted DSIFT descriptor. Dimension of words is 1024 and max pooling technique is used with pooling grid of $1\times 1,2\times 2$, and $4\times 4$. With this setup, the dimension of ScSPM feature is 21,504; this is followed by dimension reduction to $d=3000$ using PCA.

The COIL-100 dataset contains various views of 100 objects with different lighting conditions. Each object has 72 images captured from equally spaced views. Similar to the work in [80], we randomly choose 10 views of each object for training, the rest is used for test. To obtain the feature vector of each image, we first convert it to grayscale, resize to 32 $\times$ 32 pixel, vectorize this matrix to a 1024-dimensional vector, and finally normalize it to have unit norm.

Samples of this dataset are shown in Figure 18e).

Validation of efficient algorithms

To evaluate the improvement of three efficient algorithms proposed in section Discriminative dictionary learning framework, we apply these efficient algorithms and their original versions on training samples from the AR face dataset to verify the convergence speed of those algorithms. In this example, number of classes $C=100$, the random-face feature dimension $d=300$, number of training samples per class $n_{c}=n=7$, number of atoms in each particular dictionary $k_{c}=7$.

E-FDDL-D and E-FDDL-X

Figure 19 shows the cost functions and running time after each of 100 iterations of 4 different versions of FDDL: the original FDDL (O-FDDL), combination of O-FDDL-X and E-FDDL-D, combination of E-FDDL-X and O-FDDL-D, and the efficient FDDL (E-FDDL). The first observation is that O-FDDL converges quickly to a suboptimal solution, which is far from the best cost obtained by E-FDDL. In addition, while O-FDDL requires more than 12,000 seconds (around 3 hours and 20 minutes) to run 100 iterations, it takes E-FDDL only half an hour to do the same task.

E-DLSI-D and E-COPAR-D

Figure 20 and 21 compare convergence rates of DLSI and COPAR algorithms. As we can see, while the cost function value improves slightly, the run time of efficient algorithms reduces significantly. Based on benefits in both cost function value and computation, in the rest of this chapter, we use efficient optimization algorithms instead of original versions for obtaining classification results.

Visualization of learned shared bases

To demonstrate the behavior of dictionary learning methods on a dataset in the presence of shared features, we create a toy example in Figure 22. This is a classification problem with 4 classes whose basic class-specific elements and shared elements are visualized in Figure 22a). Each basis element has dimension $20$ pixel$\times 20$ pixel. From these elements, we generate 1000 samples per class by linearly combining class-specific elements and shared elements followed by noise added; 200 samples per class are used for training, 800 remaining images are used for testing. Samples of each class are shown in Figure 22b).

Figure 22c) show sample learned bases using DLSI [43] where shared features are still hidden in class-specific bases. In LC-KSVD bases (Figure 22d) and e)), shared features (the squared in the middle of a patch) are found but they are classified as bases of class 1 or class 2, diminishing classification accuracy since most of test samples are classified as class 1 or 2. The same phenomenon happens in FDDL bases (Figure 22f)).

The best classification results happen in three shared dictionary learnings (COPAR [3] in Figure 22g), JDL [33] in Figure 22h) and the proposed LRSDL in Figure 22i)) where the shared bases are extracted and gathered in the shared dictionary. However, in COPAR and JDL, shared features still appear in class-specific dictionaries and the shared dictionary also includes class-specific features. In LRSDL, class-specific elements and shared elements are nearly perfectly decomposed into appropriate sub dictionaries. The reason behind this phenomenon is the low-rank constraint on the shared dictionary of LRSDL. Thanks to this constraint, LRSDL produces perfect results on this simulated data.

Effect of the shared dictionary sizes on overall accuracy

We perform an experiment to study the effect of the shared dictionary size on the overall classification results of three shared dictionary methods: COPAR [3], JDL [33] and LRSDL in the AR gender dataset. In this experiment, 40 images of each class are used for training. The number of shared dictionary bases varies from 10 to 80. In LRSDL, because there is a regularization parameter $\eta$ which is attached to the low-rank term (see equation (33)), we further consider three values of $\eta$: $\eta=0$, i.e. no low-rank constraint, $\eta=0.01$ and $\eta=0.1$ for two different degrees of emphasis. Results are shown in Figure 23.

We observe that the performance of COPAR heavily depends on the choice of $k_{0}$ and its results worsen as the size of the shared dictionary increases. The reason is that when $k_{0}$ is large, COPAR tends to absorb class-specific features into the shared dictionary. This trend is not associated with LRSDL even when the low-rank constraint is ignored ($\eta=0$), because LRSDL has another constraint ($\|\mathbf{X}^{0}-\mathbf{M}^{0}\|_{F}^{2}$ small) which forces the coefficients corresponding to the shared dictionary to be similar. Additionally, when we increase $\eta$, the overall classification of LRSDL also gets better. These observations confirm that our two proposed constraints on the shared dictionary are important, and the LRSDL exhibits robustness to parameter choices. For JDL, we also observe that its performance is robust to the shared dictionary size, but the results are not as good as those of LRSDL.

Overall Classification Accuracy

Table 8 shows overall classification results of various methods on all presented datasets in terms of mean \pm\leavevmode\nobreak\standard deviation. It is evident that in most cases, three dictionary learning methods with shared features (COPAR [3], JDL [33] and our proposed LRSDL) outperform others with all five highest values presenting in our proposed LRSDL. Note that JDL method represents the query sample class by class. We also extend this method by representing the query sample on the whole dictionary and use the residual for classification as in SRC. This extended version of JDL is called JDL*.

Performance vs. size of training set

Real-world classification tasks often have to contend with lack of availability of large training sets. To understand training dependence of the various techniques, we present a comparison of overall classification accuracy as a function of the training set size of the different methods. In Figure 25, overall classification accuracies are reported for first five datasets888For the COIL-100 dataset, number of training images per class is already small (10), we do not include its results here. corresponding to various scenarios. It is readily apparent that LRSDL exhibits the most graceful decline as training is reduced. In addition, LRSDL also shows high performance even with low training on AR datasets.

Performance of LRSDL with varied parameters

Figure 24 shows the performance of LRSDL on the AR face dataset with different values of $\lambda_{1}$, $\lambda_{2}$, and $\eta$ with other parameters fixed. We first set these three parameters to $0.003$ then vary each parameter from $10^{-4}$ to $0.3$ while two others are fixed. We observe that the performance is robust to different values with the accuracies being greater than 98% in most cases. It also shows that LRSDL achieves the best performance when $\lambda_{1}=0.01,\lambda_{2}=0.003,\eta=0.003$.

Run time of different dictionary learning methods

Finally, we compare training and test time per sample of different dictionary learning methods on the Oxford Flower dataset. Note that, we use the efficient FDDL, DLSI, COPAR in this experiment. Results are shown in Table 9. This result is consistent with the complexity analysis reported in Table 7 with training time of LRSDL being around half an hour, 10 times faster than COPAR [3] and also better than other low-rank models, i.e. $D^{2}L^{2}R^{2}$ [45], and SRRS [80].

Discussion and Conclusion

In this chapter, our primary contribution is the development of a discriminative dictionary learning framework via the introduction of a shared dictionary with two crucial constraints. First, the shared dictionary is constrained to be low-rank. Second, the sparse coefficients corresponding to the shared dictionary obey a similarity constraint. In conjunction with discriminative model as proposed in [38, 4], this leads to a more flexible model where shared features are excluded before doing classification. An important benefit of this model is the robustness of the framework to size ($k_{0}$) and the regularization parameter ($\eta$) of the shared dictionary. In comparison with state-of-the-art algorithms developed specifically for these tasks, our LRSDL approach offers better classification performance on average.

In Section Efficient solutions for optimization problems and Efficient solutions for other dictionary learning methods, we discuss the efficient algorithms for FDDL [4], DLSI [43], then flexibly apply them into more sophisticated models. Thereafter in Section Complexity analysis and  Validation of efficient algorithms, we both theoretically and practically show that the proposed algorithms indeed significantly improve cost functions and run time speeds of different dictionary learning algorithms. The complexity analysis also shows that the proposed LRSDL requires less computation than competing models.

As proposed, the LRSDL model learns a dictionary shared by every class. In some practical problems, a feature may belong to more than one but not all classes. Very recently, researchers have begun to address this issue [84, 85]. In future work, we will investigate the design of hierarchical models for extracting common features among classes.

Classifying Multi-channel UWB SAR Imagery via Tensor

Sparsity Learning Techniques

Introduction

Over the past two decades, the U.S. Army has been investigating the capability of low-frequency, ultra-wideband (UWB) synthetic aperture radar (SAR) systems for the detection of buried and obscured targets in various applications, such as foliage penetration [86], ground penetration [87], and sensing-through-the-wall [88]. These systems must operate in the low-frequency spectrum spanning from UHF frequency band to L band to achieve both resolution and penetration capability. Although a lot of progress has been made over the years, one critical challenge that low-frequency UWB SAR technology still faces is discrimination of targets of interest from other natural and manmade clutter objects in the scene. The key issue is that that the targets of interest are typically small compared to the wavelength of the radar signals in this frequency band and have very low radar cross sections (RCSs). Thus, it is very difficult to discriminate targets and clutter objects using low-resolution SAR imagery.

SAR geometry and image formation overview

Figure 26 shows a typical side-looking SAR geometry where a vehicle or an airborne-based radar transmits wideband signals to the imaging area and measures the backscatter signals as it moves along an aperture. Let $s_{k}$ be the received range-compressed data received at $k^{\text{th}}$ aperture; the backprojection algorithm [89] computes the value for each SAR image pixel from the imaging area as follows:

[TABLE]

where $N_{i_{2}}-N_{i_{1}}+1$ is the total number of aperture records used for the coherent integration, $w_{k}$ is the weighting value at the $k^{\text{th}}$ aperture position, and $f(i,k)$ is the shift index to the signal $s_{k}$ for $i^{\text{th}}$ pixel at the $k^{\text{th}}$ aperture position. For a typical SAR image with $0^{\circ}$ aspect angle, $N_{i_{1}}$ and $N_{i_{2}}$ correspond to the angle values of $-\alpha/2$ and $\alpha/2$ to form a SAR image with an integration angle of $\alpha$. Note that for a true constant integration angle SAR image formation, $N_{i_{1}}$ and $N_{i_{2}}$ are computed for every pixel of the SAR image. However, for computational efficiency, a large image area is divided into smaller subimages. For each subimage, SAR image formation is computed using the $N_{i_{1}}$ and $N_{i_{2}}$ values derived from the geometry of the center of the subimage and the radar aperture. To exploit the aspect dependence information of target, instead of forming a single SAR image with $0^{\circ}$ aspect angle as described above with a single sector $[-\alpha/2,\alpha/2]$, we form multiple sectors to generate the corresponding SAR images at different aspect angles. For example, in our consecutive multi-look experiment, three SAR images are formed: left side image that covers sector $[-\alpha/2,\alpha/6]$, broadside image that covers sector $[-\alpha/6,\alpha/6]$, and right side image at sector $[\alpha/6,\alpha/2]$.

To achieve a constant cross-range resolution SAR image, a large image area is also divided into smaller subimages. Each subimage is formed using a different subset of aperture. For a constant integration angle, subimages at farther range would be integrated using a longer aperture than near-range subimages. The subimages are then mosaicked together to form a single large image that covers the area of interest. In consecutive multi-look SAR image processing mode, instead of generating a single SAR image, multiple SAR images are formed at different viewing angles to exploit the aspect angle dependent information from targets. Thus, each aperture for each subimage is further divided into smaller segments (either disjoint or overlapped) that are used to form consecutive multi-look SAR images as illustrated in Figure 26. The aspect angle dependence features from targets have been exploited in past research before using different techniques [90, 91, 89].

Closely related works and motivation

UWB radar techniques have recently attracted increasing attention in the area of penetration and object detection, thanks to their usage in security applications and surveillance systems [92]. T. Sakamoto et al. [93] proposed fast methods for ultra-wideband (UWB) radar imaging that can be applied to a moving target. The technology has been also applied to 3-D imaging applications [94], human posture [95], human activity [96], vital sign [97], and liquid material [98] classification problems. In these papers, due to high dimensionality and small signal-to-ratio (SNR), the signals need to be preprocessed, e.g., dimensionality reduction and background subtraction, before being used to train a classifier. It has been shown that support vector machines usually provide the best performance [96, 98]. It is worth noting that in the aforementioned applications, objects are usually big (human) and captured from a relatively small distance. On the contrary, objects in our problem are relatively small and quite far from the radar.

In this chapter, we consider the problem of discriminating and classifying buried targets of interest (metal and plastic mines, 155-mm unexploded ordinance [UXO], etc.) from other natural and manmade clutter objects (a soda can, rocks, etc.) in the presence of noisy responses from the rough ground surfaces for low-frequency UWB 2-D SAR images. For classification problems, sparse representation-based classification [1] (SRC) has been successfully demonstrated in other imagery domains such as medical image classification [19, 47, 2, 20], hyperspectral image classification [21, 22, 23], high-resolution X-band SAR image classification [24], video anomaly detection [26], and several others [73, 99, 5, 100, 27, 28, 29, 30]. However, in the low-frequency RF UWB SAR domain, although we have studied the feasibility of using SRC for higher-resolution 3-D down-looking SAR imagery [101], the application of SRC to low-frequency UWB 2-D SAR imagery has not been studied to date due to the aforementioned low-resolution issue. In this chapter, we generalize the traditional SRC to address target classification using either a single channel (radar polarization) or multiple channels of SAR imagery. Additionally, we further propose a novel discriminative tensor sparsity framework for multi-look multi-channel classification problem, which is naturally suitable for our problem. In sparse representations, many signals can be expressed by a linear combination of a few basic elements taken from a “dictionary”. Based on this theory, SRC [1] was originally developed for robust face recognition. The main idea in SRC is to represent a test sample as a linear combination of samples from the available training set. Sparsity manifests because most of the nonzero components correspond to basic elements with the same class as the test sample.

Multi-channel SRC has been investigated before in medical images [47, 20]. In these papers, one dictionary for each channel is formed from training data with locations of all channels of one training point being the same in all dictionaries. Then intuitively, when sparsely encoding each channel of a new test point using these dictionaries, we obtain sparse codes whose active (nonzero) elements tend to happen at the same locations in all channels. In other words, active elements are simultaneously located at the same location across all channels. This intuition is formulated based on $l_{0}$ pseudo-norm, which is solved using a modified version of simultaneous orthogonal matching pursuit (SOMP) [102]. The cost function is nonconvex, and hence, it is difficult to find the global solution. Furthermore, when more constraints involved, there is no straightforward way to extend the algorithm. In this chapter, we proposed another way of formulating the simultaneity constraint based on the $l_{12}$ norm, which enforces the row sparsity of the code matrix (in tensor form, we call it tube sparsity). The newly convex optimization problem can be solved effectively using the fast iterative shrinkage thresholding algorithm (FISTA) [42]. We also propose other different tensor sparsity models for our multi-channel classification problems.

It has been shown that learning a dictionary from the training samples instead of concatenating all of them as a dictionary can further enhance performance of sparsity-based methods. On one hand, the training set can be compacted into a smaller dictionary, reducing computational burden at the test time. On the other hand, by using dictionary learning, discriminative information of different classes can be trained via structured discriminative constraints on the dictionary as well as the sparse tensor code. A comprehensive study of discriminative dictionary learning methods with implementations is presented at [5, 103]. These dictionary learning methods, however, are all applied to single-channel problem where samples are often represented in form of vectors. While multi-channel signals can be converted to a long vector by concatenating all channels, this trivial modification not only leads to the curse of dimensionality of high-dimensional space, but also possibly neglects cross-channel information, which might be crucial for classification. In this chapter, we also propose a method named TensorDL, which is a natural extension of single-channel dictionary learning frameworks to multi-channel dictionary learning ones. Particularly, the cross-channel information will be captured using the aforementioned simultaneity constraint.

Naturally, when a radar carried by a vehicle or aircraft moves around an object of interest, it can capture multiple consecutive views of that object (see Fig. 26). Consequently, if the multi-look information is exploited, the classification accuracy will be improved. While the multi-look classification problem has been approached before by SRC-related methods [24, 99], none of these works uses the relative continuity of different views. We propose a framework to intuitively exploit this important information. More importantly, the optimization problem corresponding to this structure can be converted to the simultaneous sparsity model by using an elegant trick that we call ShiftSRC. Essentially, a tensor dictionary is built by circularly shifting an appropriate amount of a single-channel dictionary. When we sparsely code the multi-look signals using this tensor dictionary, the tensor sparse code becomes tube sparsity.

Contributions

The main contributions of this chapter are as follows:

A framework for simultaneously denoising and classifying 2-D UWB SAR imagery999The preliminary version of this work was presented in IEEE Radar Conference, 2017 [104]. Subtle features from targets of interest are directly learned from their SAR imagery. The classification also exploits polarization diversity and consecutive aspect angle dependence information of targets. 2. 2.

A generalized tensor discriminative dictionary learning (TensorDL) is also proposed when more training data involved. These dictionary learning frameworks are shown to be robust even with high levels of noise. 3. 3.

A relative SRC framework (ShiftSRC) is proposed to deal with multi-look data. Low-frequency UWB SAR signals are often captured at different views of objects, depending on the movement of the radar carriers. These signals contain uniquely important information of consecutive views. With ShiftSRC, this information will be comprehensively exploited. Importantly, a solution to the ShiftSRC framework can be obtained by an elegant modification on the training dictionary, resulting in a tensor sparse coding problem, which is similar to a problem proposed in contribution 1).

The remainder of this chapter is organized as follows. Section II presents different tensor sparsity frameworks and the discriminative tensor dictionary learning scheme for multi-channel classification problems. The ShiftSRC for multiple-relative-look and solutions to all proposed frameworks are also presented in this section. Section III shows extensive experimental results on a simulated dataset for several scenarios. An experiment with a realistic dataset is also included. Section IV concludes the chapter.

Sparse representation-based classification

Notation

Scalars are denoted by italic letters and may be either lower or uppercase, e.g., $d,N,k$. Vectors and matrices are denoted by bold lowercase ($\mathbf{x,y}$) and bold upper case ($\mathbf{X,Y}$), respectively. In this chapter, we also consider 3-D tensors (tensors for short) whose dimensions are named row, column, and channel. A tensor with only one column will be denoted by a bold, lowercase, calligraphic letter ($\mathbcal{x},\mathbcal{y}$). Tensors with more than one column will be denoted by an bold, uppercase, calligraphic letters ($\mathbfcal{X},\mathbfcal{Y},\mathbfcal{D}$).

For any tensor $\mathbfcal{M}$, let $\mathbfcal{M}^{(t)}$ be its $t$-th channel. For convenience, given two tensors $\mathbfcal{M},\mathbfcal{N}$, the tensor multiplication $\mathbfcal{P}=\mathbfcal{M}\mathbfcal{N}$ is considered channel-wise multiplication, i.e., $\mathbfcal{P}^{(t)}=\mathbfcal{M}^{(t)}\mathbfcal{N}^{(t)}$. For a tensor $\mathbfcal{M}$, we also denote the sum of square of all elements by $\|\mathbfcal{M}\|_{F}^{2}$ and the sum of absolute values of all elements by $\|\mathbfcal{M}\|_{1}$. Tensor addition/subtraction simply represents element-wise addition/subtraction. Each target sample is represented by a UWB SAR image formed using either a single (using co-pol) or multiple polarization (using both co-pol and cross-pol) channels. Thus, one target sample is denoted by $\mathbcal{y}\in\mathbb{R}^{d\times 1\times T}$, where $d$ is the total number of image pixels and $T$ is the number of polarization channels. A collection of $N$ samples is denoted by $\mathbfcal{Y}\in\mathbb{R}^{d\times N\times T}$.

Consider a general classification problem with $C$ different classes. Let $\mathbfcal{D}_{c}(1\leq c\leq C)$ be the collection of all training samples from class $c$, $\mathbfcal{D}_{0}$ be the collection of samples in the shared class, and $\mathbfcal{D}=[\mathbfcal{D}_{1},\dots,\mathbfcal{D}_{C},\mathbfcal{D}_{0}]$ be the total dictionary with the concatenation being done at the second dimension (column). In our problem, the shared class can be seen as the collection of ground images.

Classification scheme

Using the definition of tensor multiplication, a sparse representation of $\mathbcal{y}$ using $\mathbfcal{D}$ can be obtained by solving

[TABLE]

where $\lambda$ is a positive regularization parameter and $g(\mathbcal{x})$ is a function that encourages $\mathbcal{x}$ to be sparse. Denote by $\mathbcal{x}^{i}$ the sparse coefficient of $\mathbcal{y}$ on $\mathbfcal{D}_{i}$. Then, the tensor $\mathbcal{x}$ can be divided into $C+1$ tensor parts $\mathbcal{x}^{1},\mathbcal{x}^{2},\dots,\mathbcal{x}^{C},\mathbcal{x}^{0}$.

After solving the sparse coding problem (79), shared features (grounds in our problem) are eliminated by taking $\bar{\mathbcal{y}}=\mathbcal{y}-\mathbfcal{D}_{0}\mathbcal{x}^{0}$. Then the identity of one sample $\mathbcal{y}$ can be determined by the dictionary that provides the minimum residual:

[TABLE]

Confuser detection: In practical problems, the set of confusers is not limited to the training set. A confuser can be anything that is not a target; it can be solely the ground or a capture of an unseen object. In the former case, the test signal can be well represented by using only the ground $\mathbfcal{D}_{0}$, while in the latter case, the sparse representation assumption is no longer valid. Therefore, one test signal $\mathbcal{y}$ is classified as a confuser if one of following three conditions is satisfied: i) it is not sparsely interpreted by the total dictionary $\mathbfcal{D}$; ii) it has the most active elements in the sparse code locating at $\mathbcal{x}^{0}$; and iii) it is similar to known confusers.

Generalized sparse representation-based classification

In SRC [1] where only one channel is considered, $g(\mathbcal{x})$ is simply a function that forces sparsity as $l_{0}$- or $l_{1}$-minimization. $l_{1}$-minimization is often used since it leads to a convex optimization problem with tractable solutions. In the following, we present two natural extensions of SRC for multi-channel cases. We also proposed two tensor sparse representation methods that can enhance classification accuracy by exploiting cross-channel information. These four generalized SRC methods are as follows.

a) Apply a sparsity constraint on each column of the coefficient tensor $\mathbcal{x}$, no cross-channel constraint.

We can see that this is similar to solving $T$ separate sparse coding problems, each for a polarization channel. The classification rule is executed based on the sum of all the squares of the residuals. We refer to this framework as SRC-cumulative residual or SRC-CR (CR in the short version). (See Figure 27a with sparse tensor $\mathbcal{x}_{1}$). The sparse code corresponding to the $t^{\text{th}}$ channel, $\mathbcal{x}^{(t)}$, is computed via the traditional $l_{1}$-minimization:

[TABLE]

which can be solved effectively using FISTA [42], ADMM [44], etc., algorithms or the SPAMS toolbox [67].

b) Concatenate all channels.

The most convenient way to convert a multi-channel problem to a single-channel problem is to concatenate all $T$ channels of one signal to obtain a long vector. By doing so, we have the original SRC framework with $l_{1}$-norm minimization. After solving the sparse coding problem, if we break the long vectors and rearrange them back to tensor forms, then the tensor sparse coefficients $\mathbcal{x}$ can be formed by replicating the one-channel sparse coefficient at all channel (see Figure 27b with sparse tensor $\mathbcal{x}_{2}$). From this tensor viewpoint, the code tensor $\mathbcal{x}$ will have few active “tubes”; moreover, all elements in a tube are the same. We refer to this framework as SRC-concatenation or SRC-CC (CC in the short version).

The optimization problem in SRC-CC and its solution are very straightforward. First, we stack all channel dictionaries into a long one: $\hat{\mathbfcal{D}}=[\mathbfcal{D}^{(1)};\mathbfcal{D}^{(2)};\dots;\mathbfcal{D}^{(T)}]$ (symbol ‘;’ represents the concatenation in the first dimension). Then for every test signal, we also stack all of its channels to form a long vector: $\hat{\mathbcal{y}}=[\mathbcal{y}^{(1)};\mathbcal{y}^{(2)};\dots;\mathbcal{y}^{(T)}]$. The optimization problem (79) becomes the traditional $l_{1}$ regularization problem:

[TABLE]

then also can be solved by aforementioned methods.

c) Use a simultaneous model.

Similar to the SHIRC model proposed in [47, 20], we can impose one constraint on active elements of tensor $\mathbcal{x}$ as follows: $\mathbcal{x}$ also has few nonzero tubes as in SRC-CC; however, elements in one active tube are not necessarily the same. In other words, the locations of nonzero coefficients of training samples in the linear combination exhibit a one-to-one correspondence across channels. If the $j$-th training sample in $\mathbfcal{D}^{(1)}$ has a nonzero contribution to $\mathbcal{y}^{(1)}$, then for $t\in\{2,\dots,T\}$, $\mathbcal{y}^{(t)}$ also has a nonzero contribution from the $j$-th training sample in $\mathbfcal{D}^{(t)}$. We refer to this framework as SRC-Simultaneous or SRC-SM (SM in the short version). (See Figure 27c with sparse tensor $\mathbcal{x}_{3}$). To achieve this requirement, we can impose on the tensor $\mathbcal{x}_{3}$ (with one column and T channels) the $l_{1,2}$-minimization constraint, which is similar to the row-sparsity constraint applied on matrices in [24].

Remarks: While SHIRC uses $l_{0}$-minimization on $\mathbcal{x}$ and applies the modified SOMP [102], our proposed SRC-SM exploits the flexibility of $l_{1,2}$-regularizer since it is convex, and easily modified when more constraints are present (e.g., non-negativity). In addition, it is more efficient especially when dealing with problems of multiple samples at input.

Concretely, the optimization problem of SRC-SM can be written in the form

[TABLE]

where $\mathbcal{x}_{k::}$ denotes the $k^{\text{th}}$ tube of the tensor code $\mathbcal{x}$ and $K$ is the total column of the dictionary $\mathbfcal{D}$, and $\text{vec}(\mathbcal{x}_{k::})$ is the vectorization of the tube $\mathbcal{x}_{k::}$. This problem is similar to the joint sparse representation (JSRC) problem proposed in [24] except that SRC-SM enforces tube-sparsity instead of the row-sparsity. Details of the algorithm that solves  are described in Section Solution to optimization problems.

d) Use a group tensor model.

Intuitively, since one object is ideally represented by a linear combination of the corresponding dictionary and the shared dictionary, it is highly likely that number of active (tensor) parts in $\mathbcal{x}$ is small, i.e., most of $\mathbcal{x}^{1},\dots,\mathbcal{x}^{C},\mathbcal{x}^{0}$ are zero tensors. This suggests us a group tensor sparsity framework as an extension of [62] that can improve the classification performance, which is referred to as SRC-GT (GT in the short version). The visualization of this idea is shown in Figure 27d.

The optimization problem of SRC-GT is similar to (83) with a slight difference in the grouping coefficients:

[TABLE]

where $C+1$ is the total number of classes (including the shared ground class), and $\text{vec}(\mathbcal{x}^{c})$ is the vectorization of the group tensor $\mathbcal{x}^{c}$. Solution to this problem will be discussed next.

The overall algorithm of generalized SRC applied to multi-channel signals in the presence of a shared class is shown in Algorithm 4.

Dictionary learning for tensor sparsity

As a natural extension, we can extend the tensor sparsity models to dictionary learning ones. Most of dictionary learning methods focus on a single-channel signal, which is not suitable for models with cross-channel information. In this work, we extend single-channel dictionary learning methods to multi-channel dictionary ones by applying aforementioned tensor sparsity constraints.

In single-channel, most of discrimination dictionary learning methods, such as FDDL [4], DLSI [43], DFDL [2], LRSDL [5], etc., have a cost function that is of the form

[TABLE]

where $\bar{g}(\mathbf{X})$ is a function enforcing the sparsity of $\mathbf{X}$, and $\bar{f}_{\mathbf{Y}}(\mathbf{D},\mathbf{X})$, which includes fidelity and discriminant terms, is a function of $\mathbf{D},\mathbf{X}$ and depends on the training samples $\mathbf{Y}$.

One straightforward extension of these single-channel models to a multi-channel case is to apply the first term $f_{\mathbf{Y}}(\mathbf{D},\mathbf{X})$ to each channel and join all channels by a sparsity constraint represented by $g(\mathbf{X})$. Naturally, $g(\mathbf{X})$ can be one of four presented cross-channel sparsity constraints. Concretely, the overall cost function would be in the form

[TABLE]

where $f_{\mathbfcal{Y}^{(t)}}(\mathbfcal{D}^{(t)},\mathbfcal{X}^{(t)})=\bar{f}_{\mathbfcal{Y}^{(t)}}(\mathbfcal{D}^{(t)},\mathbfcal{X}^{(t)})$ and $g(\mathbfcal{X})$ is one of {CR, CC, SM, GT} sparsity constraints.

In this chapter, we particularly focus on extending FDDL [4] to the multi-channel case. FDDL is a special case of LRSDL [5] without an explicit shared dictionary. FDDL is chosen rather than LRSDL since in our problem, the shared dictionary, e.g., grounds, is already separated out. We also adopt fast and efficient algorithms proposed in the dictionary learning toolbox DICTOL [103] to update each channel of the dictionary $\mathbfcal{D}$. Also, the sparse code tensor $\mathbfcal{X}$ is updated using FISTA [42] algorithm, which is discussed in Section Solution to optimization problems.

The proposed cross-channel dictionary learning method is named TensorDL suffixed by CR, CC, SM, or GT when different sparsity constraints are applied on $\mathbfcal{X}$.

Tensor sparsity with multiple relative looks

In realistic situations, objects might be captured at different angles instead of only one. Moreover, these are often consecutive angles that a plane or a vehicle can capture (see Figure 28a left) while moving around objects. Based on this observation, in the training phase, we collect data from different views of objects, as in Figure 28a right, and orderly arrange them into dictionaries for each object, as illustrated in Figure 28b. The test objects, which are often captured at different relative views $\mathbf{y}_{1},\mathbf{y}_{2},\mathbf{y}_{3}$, are then sparsely coded by the whole dictionary.

Intuitively, if $\mathbf{x}_{1}$ is the sparse code of the first view $\mathbf{y}_{1}$ with only few active elements, then the sparse code $\mathbf{x}_{2}$ of the next view $\mathbf{y}_{2}$ will be active at locations shifted by one. Similarly, active locations in $\mathbf{x}_{3}$ of the third view will be shifted by two compared to $\mathbf{x}_{1}$, and so on. In this case, active coefficients form a “stair”, as illustrated in Figure 28b right.

The sparse coding problem with the “stair” sparsity is a novel problem and has not been addressed before. In this chapter, we propose a method called ShiftSRC to convert this problem to a previously known problem. The central idea is that if we stack all views into a tensor and “circularly shift” the ordered dictionary by one to the left at each view, then the obtained tensor code will be a tensor with active elements forming a “tube” sparsity (see Figure 28c). The solution to this problem is similar to the solution of SRC-SM as stated in this chapter.

Solution to optimization problems

Both optimization problems (83) and (84) have the form

[TABLE]

where

• •

$g(\mathbcal{x})$ is sum of norm 2, then it is a continuous convex function and nonsmooth.

• •

$f(\mathbcal{x})=\frac{1}{2}\|\mathbcal{y}-\mathbfcal{D}\mathbcal{x}\|_{F}^{2}$ is a continuous convex function of the type $C^{1,1}$, i.e., continuously differentiable with a Lipschitz continuous gradient $L$:

[TABLE]

We observe that, with these properties, the optimization problem (88) can be solved by FISTA [42]. There are three main tasks in FISTA:

1. Calculating $\nabla f(\mathbcal{x})$, which can be easily computed as

[TABLE]

where each channel of $\mathbfcal{D}^{T}$ is the transpose of the corresponding channel of $\mathbfcal{D}$.

2. Calculating a Lipschitz constant of $\nabla f(\mathbcal{x})$. For our function $f$, we can choose

[TABLE]

where $\lambda_{\max}$ is the maximum eigenvalue of a square matrix.

3. Solving a suboptimization problem of the form

[TABLE]

with $\eta=\frac{\lambda}{L}$.

For SRC-SM, problem (90) has the form

[TABLE]

Each problem in (91) is a minimum $l_{2}$-norm minimization with solution being

[TABLE]

Similarly, for SRC-GT, problem (90) can be written as

[TABLE]

with solution being:

[TABLE]

A step by step description of SRC-SM and SRC-GT algorithms are given in Algorithm 2.

Experimental results

In this section, we apply the proposed models to the problem of classifying objects of interest. Extensive results are presented on simulated and real-life datasets. A MATLAB toolbox for the tensor sparsity methods presented in this chapter is available at [105].

Electromagnetic (EM) Simulation data

SRC is applied to a SAR database consisting of targets (metal and plastic mines, 155-mm unexploded ordinance [UXO], etc.) and clutter objects (a soda can, rocks, etc.) buried under rough ground surfaces. The electromagnetic (EM) radar data is simulated based on the full-wave computational EM method known as the finite-difference, time-domain (FDTD) software [106], which was developed by the U.S. Army Research Laboratory (ARL). The software was validated for a wide variety of radar signature calculation scenarios [107, 108]. Our volumetric rough ground surface grid – with the embedded buried targets – was generated by using the surface root- mean-square (rms) height and the correlation length parameters. The targets are flush buried at 2-3 cm depth. In our experiments, the easiest case of rough ground surface in our experiments, the surface rms is 5.6 mm and the correlation length is 15 cm. The SAR images of various targets and clutter objects are generated from EM data by coherently integrated individual radar return signals along over a range of aspect angles. The SAR images are formed using the backprojection image formation [109] with an integration angle of $30^{\circ}$. Figure 29a shows the SAR images (using vertical transmitter, vertical receiver – VV – polarization) of some targets that are buried under a perfectly smooth ground surface. Each target is imaged at a random viewing aspect angle and an integration angle of $30^{\circ}$. Figures 29b and  29c show the same targets as Figure 29a, except that they are buried under a rough ground surface (the easiest case corresponds to ground scale/noise level101010Note that a higher ground scale means more noisy images. Henceforth, in the text as well as figures we simply refer to ground scale as noise level for ease of exposition. = 1 and harder case corresponds to noise level = 5). Similarly, Figures 29d, 29e, and 29f show the SAR images of some clutter objects buried under a smooth and rough surface, respectively. For training, the target and clutter object are buried under a smooth surface to generate high signal-to-clutter ratio images. We include 12 SAR images that correspond to 12 different aspect angles ($0^{\circ},30^{\circ},\dots,330^{\circ}$) for each target type. For testing, the SAR images of targets and confusers are generated at random aspect angles and buried under rough ground surfaces. Various levels of ground surface roughness are simulated by selecting different ground surface scaling factors when embedding the test targets under the rough surfaces. Thus, the resulting test images are very noisy with a very low signal-to-clutter ratio. Each image is a polarization signal of object which is formed by one of transmitter-receiver setups: vertical-vertical (VV), horizontal-horizontal (HH), or horizontal-vertical (HV). Each data sample of one object is represented by either i) one SAR image using data from one co-pol (VV, HH) channel or ii) two or more images using data from co-pol (VV, HH) and cross-pol (HV) channels. For each target type, we tested 100 image samples measured at random aspect angles.

Denoised signal visualization

We construct the dictionary $\mathbfcal{D}=\left[\begin{matrix}\mathbfcal{D}_{t},\mathbfcal{D}_{c},\mathbfcal{D}_{c}\end{matrix}\right]$ of all three polarizations, with $\mathbfcal{D}_{t},\mathbfcal{D}_{c},\mathbfcal{D}_{g}$ being dictionaries of targets, confusers, and grounds, respectively. The sub-dictionary $\mathbfcal{D}_{o}=\left[\begin{matrix}\mathbfcal{D}_{t},\mathbfcal{D}_{c}\end{matrix}\right]$ can be seen as the dictionary of the objects of interest. For a noisy signal $\mathbcal{y}$, we first solve the following problem:

[TABLE]

where $g(\mathbcal{x})$ is the proposed SM constraint. The sparse tensor code $\mathbcal{x}$ is then decomposed into two parts, $\mathbcal{x}^{o}$ and $\mathbcal{x}^{\text{g}}$. The latter can be considered coefficients corresponding to the ground dictionary $\mathbfcal{D}_{g}$. The original signal can be approximately decomposed into two parts: $\mathbfcal{D}_{g}\mathbcal{x}^{g}$ as separated noise, and $\mathbfcal{D}_{o}\mathbcal{x}^{o}$ as the denoised signal. Visualization of these signals are shown in Figure 30. We can see that the tensor framework successfully decomposes noisy signals into a clean part and a ground signal. These results show the potential of the proposed tensor sparsity frameworks for classification.

Overall classification accuracy

We apply four methods presented in Generalized sparse representation-based classification to different combinations of three polarizations: VV, HH, VV+HH, HH+HV, and VV+HH+HV (or ALL), and also compare these results with those obtained by a support vector machine (SVM) using the libsvm library [69]. SVM was also applied to classify UWB signals [98, 96]. The training set comprises all five target sets, four out of five confuser sets (each time we leave one confuser set out, which is meant to be unseen), and the ground set. While all confuser sets can be seen as one class – confuser class – there are two ways of organizing target sets. First, each of five targets is considered one class in case the identity of each target is crucial (we name this scenario separate-target with five target classes and one confuser class). Second, if we only need to know whether an object of interest is a target or not, we can consider all five targets as one class and its corresponding name is all-target with one target class and one confuser class. We also consider two families of the sparse coding problem, one with and one without the non-negativity constraint on the sparse code $\mathbcal{x}$ in each of the tensor sparse coding methods. Each experiment is conducted 10 times and their results are reported as the average number. Parameters in each method are chosen by a 5-fold cross-validation procedure. In this experiment, all test samples are corrupted by small noise, i.e., the noise level is set to one. In our experiments, we use overall classification accuracy as the performance metric, which computes percentage of correctly classified samples over total test samples across all classes.

Overall classification accuracies of different methods on different polarization combinations are reported in Table 10. From the table, a few inferences can be made:

• •

SRC-related methods with non-negative sparse coefficients perform better than those without this constraint111111Based on this observation, from now on, all other results are implicitly reported with the non-negativity constraint.. In addition, SVM is outperformed by all other methods in all tasks.

• •

SRC-SM provides the best results in all combinations for both the separate-target and all-target scenarios. The best accuracies in both cases are significantly high with slightly over 90% in the six-class classification problem and nearly 97% in the binary classification problem.

• •

If only one polarization is used, SRC-CR, SRC-CC, and SRC-SM have identical results, since all of them are basically reduced to traditional SRC in this case. Additionally, these three methods slightly outperform SRC-GT in this scenario.

• •

If only one polarization can be obtained, HH always outperforms VV and by a significant margin. Additionally, the HH+VV combination often worsens the results versus using HH alone.

• •

If the same method is used on different combinations, the best results are mostly obtained by the combination of HH and HV polarizations in both the separate-target and all-target scenarios.

Effect of noise levels on overall accuracy

The results in the previous section are collected in the presence of small corruption (noise level is only 1). In real problems, objects of interest are, deliberately or not, buried under extremely rough surfaces in order to fool classifiers. In this section, we conduct an experiment to see how each method performs when the level of corruption increases in the separate-target scenario.

Classification results of five different methods on different polarization combinations and different noise levels are shown in Figure 31. First of all, similar trends to small corruption can be observed in that SRC-SM shows best performance in all cases with bigger gaps occurring at high levels of noise. In addition, of the four SRC-related methods, SRC-CC is beaten by all three others when more than one polarization involved. This can be explained by the fact that SRC-CC suffers from the curse of dimensionality when each sample is represented by concatenating long vectors. It can also be seen that SRC-GT performs better than SRC-CR and SRC-CC in the presence of multiple polarizations. Last but not least, the best classification results can be seen at HH+HV and ALL among the five different combinations considered.

Effect of tensor dictionary learning on overall accuracy

We report the classification results for the all-target scenario with different noise levels. We also include the results of experiments with discriminative tensor dictionary learning methods. The results of HH+HV and ALL are reported, since they yield the best results, as observed in previous sections.

The results of nine different methods are depicted in Figure 32. These methods include SVM (the dotted line), the four SRC-related methods (dashed lines), and their corresponding tensor dictionary learning counterparts (prefixed by TensorDL, solid lines). We can see that except for the CC case, tensor dictionary learning methods outperform their corresponding SRC with gaps widening as the noise level increases. This observation confirms that tensor discriminative dictionary learning methods indeed provide accuracy improvements. Methods with the SM constraint once again emerge the winner in all noise levels, and the best accuracy numbers are observed using the HH+HV combination.

Multiple-relative-look classification accuracy

We describe a key real-world scenario where multi-look signals of an object of interest can be obtained. Figure 33a depicts the relative location of a radar carrier, a jet plane in this case, and an object of interest. The plane moves around the object at an angle corresponding to the blue arrow. One sample of object can be captured at this whole range, which can be seen as one of its looks. At the same time, the radar can also capture multiple looks of the object at smaller angles represented by green arrows. By continuing considering even smaller angles, more representatives of the object can be captured. These multiple views can provide complementary information of the object, highly likely resulting in better classification accuracy of the system.

For the training samples, for each object, we generate two sets of signals. Each set contains samples captured after each $15^{\circ}$, and each set has total of 24 views. Both sets start at the same relative angle but the first is captured by an integration angle of $30^{\circ}$, the angle in the second set is $15^{\circ}$. For the test samples, each object is represented by three signals: one by an integration angle of $30^{\circ}$ and two others by an integration angle of $15^{\circ}$, as depicted in Figure 33a. Similar to previous experiments, test samples are captured at random relative angles. Ground samples are also simulated in the same way. Based on three signals captured, we establish three different way of using this information in the classification process:

1look: for each object, we use signals at integration angle of $30^{\circ}$ only. If only one polarization is used, an object of interest can be identified by SRC (and implicitly, SVM). If more than one polarization is involved, one object will be determined by SRC-SM, as this is the best method based on previous experiments (see Figure 33b). 2. 2.

2look: each object is represented by two singles captured at $15^{\circ}$. This multi-look classification problem can be solved by joint SRC (JSRC) [24], or the proposed relative-look SRC (ShiftSRC) (see Figure 33c). 3. 3.

3look: uses all three signals to represent an object. In this case, the relationship between the $30^{\circ}$ signal and the first $15^{\circ}$ signal can be modeled by the SRC-SM, while the relationship between two $15^{\circ}$ signals can be formed by either JSRC or ShiftSRC (see Figure 33d).

For this experiment, we consider the all-target scenario and two polarization combinations, HH and HH+HV. It is worth noting that for the HH+HV case, there will be four channels in 2look and six channels in 3look. The results of different methods are shown in Figure 34. We can see that SVM still performs well at the lowest noise level (1), but drastically degrades with a little more corruption. On the other hand, SRC-based methods obtain good results even if the noise level is large for the HH+HV combination. Of sparse representation methods, ShiftSRC outperforms the others with the gap becoming larger for highly corrupted signals. Interestingly, ShiftSRC at 2look provides even better results than JSRC does at 3look. These results confirm the advantages of the proposed relative-look sparse representation model.

Overall accuracy on measured UWB SAR data

In this section, the results of this technique are illustrated using the data from the ARL UWB low-frequency SAR, which transmits radar signals occupying the frequency spectrum that span approximately from 50 to 1150 MHz [110]. Figure 35 shows a SAR image formed using data collected at Yuma Proving Grounds (YPG) [87]. The scene includes several rows of buried mines surrounded by clutter objects such as bushes, rocks, tire tracks, etc.

The set contains signals of 50 targets and 50 confusers. Each signal has resolution of $90\times 30$ and already includes noise from the ground. Visualization of six samples in each class are shown in Figure 36. We conduct the all-target experiment and report results of different methods on different polarization combinations in Figure 37. For each combination, three competing methods are considered: SVM, SRC, and TensorDL (both with the SM constraint). Since grounds are fixed in this data set, we report the results based on size of the training set. For each training size $N$ ($N=$10, 20, 30, or 40), we randomly choose $N$ samples from each class for training; the rest $50-N$ samples are considered test samples.

The results are reported in Figure 37 as the average of 10 tests. In general, the tensor dictionary learning performs better than sparse representation in all combinations except for the HH case. SVM also provides good results but is outperformed by other competing methods.

Conclusion

We have developed a novel discrimination and classification framework for low-frequency UWB SAR imagery using sparse representation-based methods. The framework is applicable for either single channel or multiple channel (polarizations, multilook) SAR imagery. The techniques are tested and the discrimination/classification performance of targets of interest versus natural and manmade clutter in challenging scenarios is measured using both rigorous electromagnetic simulation data and real data from the ARL UWB radar. The classification results show encouraging potential of tensor sparsity methods, even when the test images are very noisy (buried under extremely rough ground surfaces), targets have small RCSs, and the targets’ responses have very little detailed structures, i.e., targets are small compared to the wavelengths of the radar signals. The SRC-SM technique and its dictionary learning version consistently offers the best results with the combination of co- and cross-pol data, e.g., HH+HV or ALL. In addition, the non-negativity constraints on sparse tensor codes enhance the performance of the systems. Furthermore, we also show the ShiftSRC model is particularly suitable for problems with multiple-relative-look signals.

Conclusions and Future Directions

Summary of Main Contributions

The overarching theme in this dissertation is the design of signal and image classification algorithms by exploiting structurally meaningful prior information of signal representations in associated models. Different discriminative models have been explored based on signal sparsity structures. We have primarily considered low-training classification scenarios where the different signal representations exhibit discriminative structure that can be leveraged for robustness benefits.

In Chapter 2, we focus on extracting discriminative features of histopathological images that can be useful for classification and detection tasks. In particular, we based our contributions on dictionary learning to discriminative learning of complicated medical image features. By simply build one dictionary for each image class that promote small intra-class representation and prevent inter-class representation, we obtain class-specific dictionaries that explicitly capture crucial features. The framework is theoretically and experimentally shown to have a low complexity compared to other related methods. In order to verify that the framework is indeed applicable in a variety of scenarios, we conduct several experiments on three diverse histopathological datasets. It is illustrated our method is competitive with or outperforms state of the art alternatives, particularly in the regime of realistic or limited training set size. It is also shown that with minimal parameter tuning and algorithmic changes, the proposed method can be easily applied on different problems with different natures which makes it a good candidate for automated medical diagnosis instead of using customized and problem specific frameworks for every single diagnosis task. We also create a software toolbox available to help deploy the method widely as a diagnostic tool in existing histopathological image analysis systems. Particular problems such as grading and detecting specific regions in histopathology may be investigated using our proposed techniques.

Chapter 3, we propose a generalized dictionary learning method for different image classification problem. Our primary contribution in this chapter is the development of a discriminative dictionary learning framework via the introduction of a shared dictionary with two crucial constraints. First, the shared dictionary is constrained to be low-rank. Second, the sparse coefficients corresponding to the shared dictionary obey a similarity constraint. In conjunction with a widely used discriminative model, this leads to a more flexible model where shared features are excluded before doing classification. An important benefit of this model is the robustness of the framework to size and the regularization parameter of the shared dictionary. In comparison with state-of-the-art algorithms developed specifically for these tasks, the proposed model offers better classification performance on average. Another important contribution of this chapter is the dictionary learning toolbox – DICTOL. The toolbox implements many widely used generative and discriminative dictionary learning methods. Many classical methods are sped up via new numerical optimization innovations.

Table 11 summarizes different characteristics of widely used dictionary learning methods and SRC. The proposed LRSDL requires a smallest amount of training time compared with other frameworks using their original algorithms. In terms of classification performance, two methods focusing on extracting shared features, COPAR [3] and LRSDL, outperform other methods with slightly better but much faster results provided by LRSDL, thanks to the efficient algorithm and low-rank structure.

As an extension of discriminative sparsity models, we continue to propose different tensor sparsity structures for the problem of discriminating UWB-SAR imagery in Chapter 5. In this chapter, we present three novel sparsity-driven techniques, which not only exploit the subtle features of raw captured data but also take advantage of the polarization diversity and the aspect angle dependence information from multi-channel SAR data. First, the traditional sparse representation-based classification is generalized to exploit shared information of classes and various sparsity structures of tensor coefficients for multichannel data. Corresponding tensor dictionary learning models are consequently proposed to enhance classification accuracy. Lastly, a new tensor sparsity model is proposed to model responses from multiple consecutive looks of objects, which is a unique characteristic of the dataset we consider. Extensive experimental results on a high-fidelity electromagnetic simulated dataset and radar data collected from the U.S. Army Research Laboratory side-looking SAR demonstrate the advantages of proposed tensor sparsity models.

Potential Future Research Directions

The contributions in the previous chapters naturally point towards various directions for future research. We mention some of the possible extensions in this section.

Hierarchically shared feature learning

As proposed, the LRSDL model learns a dictionary shared by every class. In some practical problems, a feature may belong to more than one but not all classes. In future work, we will investigate the design of hierarchical models for extracting common features among classes.

Exploiting shared features using deep learning models

Deep learning has become a powerful framework for many image processing and computer vision applications such as image classification and object detection. To the best of our knowledge, there is no deep learning model that can reduce the effect of the shared features. As viable future research direction for this line of work, we propose to identify meaningful physical prior information for use in deep networks for extracting shared features and to demonstrate its benefits, especially in low training data scenarios.

Evolution of structured sparsity

Almost all dictionary learning methods enforce sparsity using $\ell_{1}$ norm or $\ell_{0}$ pseudo norm. It has been shown that Bayesian inference can effectively capture sparsity of signals by introducing new priors [111, 112]. Amongst those priors, a well-suited sparsity promoting prior is Spike and Slab prior which is widely used in Bayesian inference [113, 114, 115]. In fact, it is acknowledged that Spike and Slab prior is indeed the gold standard for inducing sparsity in Bayesian inference [115]. Naturally, we can incorporate this prior into a dictionary learning framework to further boost the classification performance. The main obstacle of using this prior is that coefficient update optimization problem is known to be a hard non-convex problem where most of existing solutions involve simplifying assumptions and/or relaxation. Fortunately, this hard problem has been addressed by two very recent works, ICR [16] and AMP [32]. Both ICR and AMP propose algorithms to directly solve the problem while AMP can solve the problem very effectively. These promising results can be extended into a Bayesian dictionary learning framework to further enhance the performance.

\Appendix

Complexity Analysis of Dictionary Learning methods

In this section, we compare the computational complexity for the proposed DFDL and competing dictionary learning methods: LC-KSVD[37], FDDL[38], and Nayak’s[48]. The complexity for each dictionary learning method is estimated as the (approximate) number of operations required by each method in learning the dictionary. For simplicity, we assume that number of training samples, number of dictionary bases in each class are the same, which means: $N_{i}=N_{j}=N,k_{i}=k_{j}=k,\forall i,j=1,2,\dots,c$, and also $L_{i}=L_{j}=L,\forall i,j=1,2,\dots,c.$ For the consitence, we have changed notations in those methods by denoting $\mathbf{Y}$ as training samples and $\mathbf{X}$ as the sparse code.

In most of dictionary learning methods, the complexity of sparse coding step, which is often a $l_{0}$ or $l_{1}$ minimization problem, dominates that of dictionary update step, which is typically solved by either block coordinate descent[36] or singular value decomposition[17]. Then, in order to compare the complexity of different dictionary learning methods, we focus on comparing the complexity of sparse coding steps in each iteration.

Complexity of the DFDL

The most expensive computation in DFDL is solving an Orthogonal Matching Pursuit (OMP [66]) problem. Given a set of samples $\mathbf{Y}\in\mathbb{R}^{d\times N}$, a dictionary $\mathbf{D}\in\mathbb{R}^{d\times k}$ and sparsity level $L$, the OMP problem is:

[TABLE]

R. Rubinstein et al.[116] reported the complexity of Batch-OMP when the dictionary is stored in memory in its entirety as: $T_{\text{b-omp}}=N(2dk+L^{2}k+3Lk+L^{3})+dk^{2}.$ Assuming an asymptotic behavior of $L\ll k\approx d\ll N$, the above expression can be simplified to:

[TABLE]

This result will also be utilized in analyzing complexity of LC-KSVD.

The sparse coding step in our DFDL consists of solving $c$ sparse coding problems: $\hat{\mathbf{X}}=\arg\min_{\|\mathbf{X}\|_{0}\leq L}\left\|\hat{\mathbf{Y}}-\mathbf{D}_{i}\hat{\mathbf{X}}_{i}\right\|_{F}^{2}.$ With $\hat{\mathbf{Y}}\in\mathbb{R}^{d\times cN},\mathbf{D}_{i}\in\mathbb{R}^{d\times k}$, each problem has complexity of $k(cN)(2d+L^{2})$. Then the total complexity of these $c$ problems is: $T_{\text{DFDL}}\approx c^{2}kN(2d+L^{2})$.

Complexity of LC-KSVD

We consider LC-KSVD1 only (LC-KSVD2 has a higher complexity) whose optimization problem is written as [37]:

[TABLE]

and it is rewritten in the K-SVD form:

[TABLE]

Since $\mathbf{Q}\in\mathbb{R}^{ck\times cN}$ and $\mathbf{A}\in\mathbb{R}^{ck\times ck}$, $\tilde{\mathbf{Y}}=\left[\begin{matrix}\mathbf{Y}\\ \sqrt{\alpha}\mathbf{Q}\end{matrix}\right]\in\mathbb{R}^{(d+ck)\times cN}$ and $\tilde{\mathbf{D}}=\left[\begin{matrix}\mathbf{D}\\ \sqrt{\alpha}\mathbf{A}\end{matrix}\right]\in\mathbb{R}^{(d+ck)\times ck}$. Omitting the computation of scalar multiplications, the complexity of (97) is:

[TABLE]

Complexity of Nayak’s

The optimization problem in Nayak’s[48] is:

[TABLE]

$\mathbf{X}$ is estimated via the gradient descent method that is an iterative method whose main computational task in each iteration is to calculate the gradient of $Q(\mathbf{X})=\left\|\mathbf{Y}-\mathbf{D}\mathbf{X}\right\|_{F}^{2}+\left\|\mathbf{X}-\mathbf{W}\mathbf{Y}\right\|_{F}^{2}$ with respect to $\mathbf{X}$. We have:

[TABLE]

where $\mathbf{D}^{T}\mathbf{D}+\mathbf{I}$, and $(\mathbf{D}^{T}-\mathbf{W})\mathbf{Y}$ could be precomputed and at each step, only $(\mathbf{D}^{T}\mathbf{D}+\mathbf{I})\mathbf{X}$ need to be recalculated after $\mathbf{X}$ is updated. With $\mathbf{D}\in\mathbb{R}^{d\times ck},\mathbf{X}\in\mathbb{R}^{ck\times cN},\mathbf{Y}\in\mathbb{R}^{d\times cN},\mathbf{W}\in\mathbb{R}^{ck\times d}$, the complexity of the sparse coding step can be estimated as:

[TABLE]

with $q$ being the average number of iterations needed for convergence. Here we have ignored matrix subtractions, additions and scalar multiplications and focused on matrix multiplications only. We have also used the approximation that complexity of $\mathbf{A}\mathbf{B}$ is $2mnp$ where $\mathbf{A}\in\mathbb{R}^{m\times n},\mathbf{B}\in\mathbb{R}^{n\times p}$. The first term in (98) is of $\mathbf{D}^{T}\mathbf{D}+\mathbf{I}$ (note that this matrix is symmetric, then it needs only half of regular operations), the second term is of $(\mathbf{D}^{T}-\mathbf{W})\mathbf{Y}$ and the last one comes from $q$ times complexity of calculating $(\mathbf{D}^{T}\mathbf{D}+\mathbf{I})\mathbf{X}$.

Complexity of FDDL

The sparse coding step in FDDL[38] requires solving $c$ class-specific problems:

[TABLE]

with $\mathbf{D}=[\mathbf{D}_{1},\dots,\mathbf{D}_{c}],\mathbf{X}_{i}^{T}=[(\mathbf{X}_{i}^{1})^{T},\dots,(\mathbf{X}_{i}^{c})^{T}]$, and $\mathbf{M}_{k}=[\mathbf{m}_{k},\dots,\mathbf{m}_{k}]\in\mathbb{R}^{ck\times N},\mathbf{M}=[\mathbf{m},\dots,\mathbf{m}]\in\mathbb{R}^{ck\times N}$ where $\mathbf{m}_{k}$ and $\mathbf{m}$ are the mean vector of $\mathbf{X}_{i}$ and $\mathbf{X}=[\mathbf{X}_{1},\dots,\mathbf{X}_{c}]$ respectively. The algorithm for solving this problem uses Iterative Projective Method[117] whose complexity depends on computing gradient of six Frobineous-involved terms in the above optimization problem at each iteration.

For the first three terms, the gradient could be computed as:

[TABLE]

where $\mathbf{D}^{T}\mathbf{D},\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }\mathbf{D}^{T}\mathbf{Y}_{i}$ could be precomputed with the total cost of $(ck)d(ck)+2(ck)dN=cdk(2N+ck)$; $\mathbf{D}_{i}^{T}\mathbf{D}_{i},\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }\mathbf{D}_{i}^{T}\mathbf{Y}_{i}$ could be extracted from $\mathbf{D}^{T}\mathbf{D},\text{\leavevmode\nobreak\ and\leavevmode\nobreak\ }\mathbf{D}^{T}\mathbf{Y}_{i}$ at no cost; at each iteration, cost of computing $(\mathbf{D}^{T}\mathbf{D})\mathbf{X}_{i}$ is $2(ck)^{2}N$, each of $(\mathbf{D}_{j}^{T}\mathbf{D}_{j})\mathbf{X}_{i}^{j}$ could be attained in the intermediate step of computing $(\mathbf{D}^{T}\mathbf{D})\mathbf{X}_{i}$. Therefore, with $q$ iterations, the computational cost of (Complexity of FDDL) is:

[TABLE]

For the last three terms, we will prove that:

[TABLE]

Indeed, let $\mathbf{E}_{m,n}$ be a all-one matrix in $\mathbb{R}^{m\times n}$, one could easily verify that:

[TABLE]

[TABLE]

Thus, (102) can be obtained by:

[TABLE]

For (103), with simple algebra, we can prove that:

[TABLE]

[TABLE]

Compared to (100), calculating (102), (103) and (104) require much less computation. As a result, the total cost of solving $\mathbf{X}_{i}$ approximately equals to (101); and the total estimated cost of sparse coding step of FDDL is estimated as $c$ times cost of each class-specific problem and approximately equals to:

[TABLE]

Final analyzed results of these four methods are reported in Table 1.

\Appendix

Proof of Lemmas in Chapter Fast Low-rank Shared Dictionary Learning for Image Classification

Proof of Lemma 1

Let $\mathbcal{w}_{c}\in\{0,1\}^{K}$ is a binary vector whose $j$-th element is one if and only if the $j$-th columns of $\mathbf{D}$ belong to $\mathbf{D}_{c}$, and $\mathbf{W}_{c}=\text{diag}(\mathbcal{w}_{c})$. We observe that $\mathbf{D}_{c}\mathbf{X}^{c}_{i}=\mathbf{D}\mathbf{W}_{c}\mathbf{X}_{i}$. We can rewrite $f_{\mathbf{Y},\mathbf{X}}(\mathbf{D})$ as:

[TABLE]

where we have defined:

[TABLE]

Let:

[TABLE]

From definition of $\mathbf{W}_{j}$, we observe that ‘left-multiplying’ a matrix by $\mathbf{W}_{j}$ forces that matrix to be zero everywhere except the $j$-th block row. Similarly, ‘right-multiplying’ a matrix by $\mathbf{W}_{j}^{T}=\mathbf{W}_{j}$ will keep its $j$-th block column only. Combining these two observations, we can obtain the result:

[TABLE]

Then:

[TABLE]

Lemma 1 has been proved. $\square$

Proof of Lemma 2

We need to prove two parts:

For the gradient of $f$, first we rewrite:

[TABLE]

Then we obtain:

[TABLE]

For the gradient of $g$, let $\mathbf{E}_{p}^{q}$ be the all-one matrix in $\mathbb{R}^{p\times q}$. We can verify that:

[TABLE]

We have:

[TABLE]

[TABLE]

Therefore we obtain:

[TABLE]

For $\mathbf{M}_{c}-\mathbf{M}$, first we write it in two ways:

[TABLE]

Then we infer:

[TABLE]

$\displaystyle\Rightarrow\frac{\partial}{\partial\mathbf{X}_{l}}\frac{1}{2}\sum_{c=1}^{C}\|\mathbf{M}_{c}-\mathbf{M}\|_{F}^{2}=\mathbf{M}_{l}-\mathbf{M}+\frac{1}{N}\sum_{c=1}^{C}(\mathbf{M}-\mathbf{M}_{c})\mathbf{E}_{n_{c}}^{n_{l}}.$

Now we prove that $\displaystyle\sum_{c=1}^{C}(\mathbf{M}-\mathbf{M}_{c})\mathbf{E}_{n_{c}}^{n_{l}}=\mathbf{0}$. Indeed,

[TABLE]

Then we obtain:

[TABLE]

Combining (106), (109) and $\frac{\partial\frac{1}{2}\|\mathbf{X}\|_{F}^{2}}{\partial\mathbf{X}}=\mathbf{X}$ , we have:

[TABLE]

Lemma 2 has been proved. $\square$

Proof of Lemma 3

When $\mathbf{Y},\mathbf{D},\mathbf{X}$ are fixed, we have:

[TABLE]

Let $\tilde{\mathbf{Y}}=\mathbf{Y}-\mathbf{D}\mathbf{X}$, $\hat{\mathbf{Y}}_{c}=\mathbf{Y}_{c}-\mathbf{D}_{c}\mathbf{X}_{c}^{c}$ and $\hat{\mathbf{Y}}=\left[\begin{matrix}\hat{\mathbf{Y}}_{1}&\hat{\mathbf{Y}}_{2}&\dots&\hat{\mathbf{Y}}_{C}\end{matrix}\right]$, we can rewrite (Proof of Lemma 3) as:

[TABLE]

We observe that:

[TABLE]

Now, by letting $\displaystyle\mathbf{V}=\frac{\tilde{\mathbf{Y}}+\hat{\mathbf{Y}}}{2}$, Lemma 3 has been proved.        $\square$