Deep CNNs Meet Global Covariance Pooling: Better Representation and Generalization

TL;DR

This paper introduces a robust global covariance pooling method for deep CNNs, enhancing their ability to capture rich feature statistics and improve generalization across various visual recognition tasks.

Contribution

It proposes the MPN-COV pooling method with a Gaussian embedding network and efficient normalization, addressing high-dimensional covariance estimation challenges in deep CNNs.

Findings

Outperforms existing methods on large-scale classification tasks

Achieves state-of-the-art accuracy in multiple visual recognition benchmarks

Provides a modular approach easily integrated into existing CNN architectures

Abstract

Compared with global average pooling in existing deep convolutional neural networks (CNNs), global covariance pooling can capture richer statistics of deep features, having potential for improving representation and generalization abilities of deep CNNs. However, integration of global covariance pooling into deep CNNs brings two challenges: (1) robust covariance estimation given deep features of high dimension and small sample size; (2) appropriate usage of geometry of covariances. To address these challenges, we propose a global Matrix Power Normalized COVariance (MPN-COV) Pooling. Our MPN-COV conforms to a robust covariance estimator, very suitable for scenario of high dimension and small sample size. It can also be regarded as Power-Euclidean metric between covariances, effectively exploiting their geometry. Furthermore, a global Gaussian embedding network is proposed to incorporate first-order statistics into MPN-COV. For fast training of MPN-COV networks, we implement an iterative matrix square root normalization, avoiding GPU unfriendly eigen-decomposition inherent in MPN-COV. Additionally, progressive 1x1 convolutions and group convolution are introduced to compress covariance representations. The proposed methods are highly modular, readily plugged into existing deep CNNs. Extensive experiments are conducted on large-scale object classification, scene categorization, fine-grained visual recognition and texture classification, showing our methods outperform the counterparts and obtain state-of-the-art performance.

Tables

Table 1. TABLE I: Time (ms) taken by matrix decomposition (single precision arithmetic) of a 256 × 256 256 256 256\times 256 covariance matrix.

Algorithm	CUDA cuSOLVER	Matlab (CPU function)	Matlab (GPU function)
EIG	21.3	1.8	9.8
SVD	52.2	4.1	11.9

Table 2. TABLE II: Top-1 error (%, 1-crop prediction) of GCP with various normalization methods with AlexNet on ImageNet.

method	MPN (ours)	M-Fro	M-${\ell }_2$	E-PN	top-1 Err.
Baseline	–	–	–	–	41.52
GCP	$\times$	$\times$	$\times$	$\times$	40.41
	$✓$	$\times$	$\times$	$\times$	38.51
	$✓$	$✓$	$\times$	$\times$	39.93
	$✓$	$\times$	$✓$	$\times$	39.62
	$✓$	$\times$	$\times$	$✓$	40.75
	$\times$	$✓$	$\times$	$\times$	39.87
	$\times$	$\times$	$✓$	$\times$	39.65
	$\times$	$\times$	$\times$	$✓$	39.89

Table 3. TABLE III: Impact of post-compensation on iSQRT-COV with ResNet-50 on ImageNet, evaluated with 1-crop prediction.

Pre-normalization	Post-compensation	Top-1 Err.	Top-5 Err.
Trace	w/o	N/A	N/A
	w/  BN [60]	23.12	6.60
	w/  Trace	22.14	6.22

Table 4. TABLE IV: Error rate (%, 1-crop prediction) and extra network parameters (Params.) of iSQRT-COV using various compact covariance representations (Repr.) on ImageNet. ResNet-50 is used as backbone model and is compared as baseline.

Method	DR	Repr.	Top-1/Top-5	Params.
GAP [4]	N/A	2K	24.7/7.8	25.56M
iSQRT -COV	$2048\to 256$	32K	22.14/6.22	+30.53M
	$2048\stackrel{512}{\to }256$	32K	21.72/5.99	+31.36M
	$2048\to 128$	8K	22.78/6.43	+6.26M
	$2048\stackrel{512}{\to }128$	8K	22.33/6.28	+7.07M
	$2048\to 64$	2K	23.73/6.99	+0.13M
	$2048\stackrel{256}{\to }64$	2K	22.98/6.61	+0.54M
	$4th$ row & $1G$	2K	22.40/6.35	+18.9M
	$4th$ row & $2G$	2K	22.84/6.60	+9.6M
	$4th$ row & 4G	2K	23.11/6.77	+4.3M

Table 5. TABLE V: Error rate (%, 1-crop prediction) and time of FP+BP (ms) per image of different covariance pooling methods with AlexNet on ImageNet. Numbers in parentheses indicate FP time. ∗ Following  [ 24 ] , improved B-CNN successively performs matrix square root and E-PN.

Method	Top-1 Err.	Top-5 Err.	Time
AlexNet [1]	41.8	19.2	1.32 (0.77)
B-CNN [16]	39.89	18.32	1.92 (0.83)
DeepO2P [15]	42.16	19.62	11.23 (7.04)
Improved B-CNN∗[24]	40.75	18.91	15.48 (13.04)
MPN-COV-Net	38.51	17.60	3.89 (2.59)
G2DeNet	38.48	17.59	9.86 (5.88)
iSQRT-COV-Net(Frob.)	38.78	17.67	2.56 (0.81)
iSQRT-COV-Net(trace)	38.45	17.52	2.55 (0.81)

Table 6. TABLE VI: Comparison of classification error (%, 1-crop prediction) with local second-order networks with ResNet as backbone models on ImageNet.

Method	Backbone model	Top-1 Err.	Top-5 Err.
ResNet-50 [4]	ResNet-50	24.7	7.8
FBN [46]		24.0	7.1
SORT [47]		23.82	6.72
SE-Net [62]		23.29	6.62
CBAM [63]		22.66	6.31
MPN-COV-Net		22.73	6.54
iSQRT-COV-Net		22.14	6.22
ResNet-101 [4]	ResNet-101	23.6	7.1
SE-Net [62]		22.38	6.07
CBAM [63]		21.51	5.69
iSQRT-COV-Net		21.21	5.68
ResNet-152 [4]	ResNet-152	23.0	6.7
DenseNet-201 [5]	DenseNet-201	22.58	6.34
iSQRT-COV-Net		20.69	5.48

Table 7. TABLE VII: Comparison of classification error (%, 10-crop prediction) on Places365 dataset. B-CNN and iSQRT-COV-Net use ResNet-50 as backbone model.

	VGG-VD16 [2]	GoogleNet [3]	ResNet-50 [4]	ResNet-152 [32]	B-CNN	iSQRT-COV-Net	iSQRT-COV + E-PN
Top-1 Error	44.76	46.37	44.82	45.26	44.24	43.68	45.34
Top-5 Error	15.09	16.12	14.71	14.92	14.27	13.73	15.13

Table 8. TABLE VIII: Comparison of classification accuracy (%) with state-of-the-art methods on fine-grained benchmarks.

	Method	Birds	Aircrafts	Cars
VGG-D16	VGG-VD16 [16]	70.4	76.6	79.8
	NetVLAD [49]	81.9	81.8	88.6
	NetFV [68]	79.9	79.0	86.2
	B-CNN [16]	84.0	83.9	90.6
	CBP [39]	84.3	84.1	91.2
	LRBP [40]	84.2	87.3	90.9
	KP [42]	86.2	86.9	92.4
	HIHCA [43]	85.3	88.3	91.7
	Improved B-CNN[24]	85.8	88.5	92.0
	SMSO [45]	85.01	N/A	N/A
	MPN-COV-Net	86.7	89.9	92.2
	G2DeNet	87.1	89.0	92.5
	iSQRT-COV-Net	87.2	90.0	92.5
ResNet-50	CBP [39]	81.6	81.6	88.6
	KP [42]	84.7	85.7	91.1
	SMSO [45]	85.77	N/A	N/A
	iSQRT-COV-Net (8K)	87.3	89.5	91.7
	iSQRT-COV-Net (32K)	88.1	90.0	92.8
iSQRT-COV-Net with ResNet-101		88.7	91.4	93.3

Table 9. TABLE IX: Comparison of classification accuracy ( % percent \% ) with state-of-the-art methods on DTD and Indoor67.

Method	Backbone Model	DTD	Indoor67
VGG-VD16 [41]	VGG-VD16	62.9$\pm$0.8	67.6
B-CNN [16]	VGG-VD16	72.9$\pm$0.8	79.0
FASON [41]	VGG-VD16	72.9$\pm$0.7	80.8
Deep-TEN [69]	ResNet-50	N/A	76.2
MFAFV-Net [52]	VGG-VD16	N/A	81.1
SMSO [45]	VGG-VD16	69.26	79.45
	ResNet-50	72.51	79.68
iSQRT-COV-Net	VGG-VD16	74.0$\pm$0.8	81.4
	ResNet-50	74.8$\pm$1.0	83.5

Table 10. TABLE X: Classification accuracy ( % percent \% ) of B-CNN and iSQRT-COV with different fine-tuning strategies on Birds (B), Aircrafts (A), Cars (C), DTD (D) and Indoor67 (I) datasets.

	Method	B	A	C	D	I
VGG-VD16	VGG-VD16	70.4	76.6	79.8	62.9$\pm$0.8	67.6
	B-CNN	84.0	83.9	90.6	72.9$\pm$0.8	79.0
	B-CNN (PreTr-GCP)	84.7	85.1	91.1	N/A	N/A
	iSQRT-COV (PreTr-GAP)	86.1	88.9	91.8	73.0$\pm$1.0	79.9
	iSQRT-COV + E-PN	87.1	90.1	92.7	73.2$\pm$0.9	81.3
	iSQRT-COV-Net	87.2	90.0	92.5	74.0$\pm$0.8	81.4
ResNet-50	ResNet-50 [42]	78.4	79.2	84.7	N/A	N/A
	B-CNN	84.4	85.4	91.4	N/A	N/A
	B-CNN (PreTr-GCP)	85.2	86.1	91.5	N/A	N/A
	iSQRT-COV (PreTr-GAP)	86.7	89.6	92.2	73.4$\pm$0.9	82.5
	iSQRT-COV + E-PN	88.2	90.4	93.2	74.6$\pm$0.9	82.9
	iSQRT-COV-Net	88.1	90.0	92.8	74.8$\pm$1.0	83.5

Table 11. TABLE XI: Classification error ( % percent \% ) of different methods on iNaturalist Challenge 2018.

Method	Description	Top-3 Err.
ResNet-152	320$\times$320 input	16.623
	392$\times$392 input	16.024
iSQRT-COV-Net	320$\times$320 input	15.038
	392$\times$392 input	14.704
Our Fusion I	3 $\times$ ResNet-152 Model	14.625
Our Fusion II	3 $\times$ iSQRT-COV-Net Model	13.409
Our Fusion III	Fusion I + Fusion II	13.068
Runner-up	Deep Learning Analytics	14.214

Table 12. TABLE XII: Comparison of different Gaussian embedding methods using VGG-VD16 on Birds-CUB200-2011 dataset.

Method	Gaussian Embedding	Accuracy $(\%)$
Nakayama et al. [59]	$𝐳={\left[\begin{array}{c}vec\left(𝚺+𝝁{𝝁}^T\right),{𝝁}^T\end{array}\right]}^T$	$83.5$
Calvo et al. [56] or Lovrić et al. [57]		$84.1$
[72, 73] (i.e., Calvo et al. [56] or Lovrić et al. [57] + Log-Euclidean [11])		$83.8$
G2DeNet (Ours)		$87.1$

Table 13. TABLE XIII: Comparison of our iSQRT-COV and GAP using various ResNets in terms of network parameters, floating point operations per second (FLOPs), inference time per image, memory cost for testing one image and classification error (%, 1-crop prediction), given input images of size 224 × \times 224.

Methods	Parameter	GFLOPs.	Inference time (ms)	Memory Cost	Top-1 Err.	Top-5 Err.
ResNet18 + GAP	11.69M	1.82	0.60	763M	30.24	10.92
ResNet18 + iSQRT-COV (32k)	44.20M	3.22	0.95	1263M	24.52	7.77
ResNet18 + iSQRT-COV (8k)	19.49M	3.10	0.85	883M	24.93	7.86
ResNet34 + GAP	21.80M	3.67	0.88	901M	26.70	8.58
ResNet34 + iSQRT-COV (32k)	54.31M	5.77	1.15	1407M	22.89	6.64
ResNet34 + iSQRT-COV (8k)	29.71M	5.56	1.10	1029M	23.20	6.89
ResNet50 + GAP	25.56M	3.86	1.29	975M	24.70	7.80
ResNet50 + iSQRT-COV (32k)	56.93M	6.31	1.50	1463M	22.14	6.22
ResNet50 + iSQRT-COV (8k)	32.32M	6.19	1.49	1085M	22.33	6.28
ResNet101 + GAP	44.55M	7.57	1.72	1233M	23.60	7.10
ResNet101 + iSQRT-COV (32k)	75.92M	10.02	1.87	1733M	21.21	5.68
ResNet101 + iSQRT-COV (8k)	51.31M	9.903	1.83	1355M	21.32	5.70
ResNet152 + GAP	60.19M	11.28	2.55	1483M	23.00	6.70

Table 14. TABLE XIV: Classification accuracy (%) of iSQRT-COV-Net with the softmax and SVM classifiers on fine-grained visual classification task.

	Backbone	Classifier	Birds	Aircrafts	Cars
iSQRT-COV -Net	VGG-D16	SVM	87.2	90.0	92.5
		softmax	87.1	90.0	92.5
	ResNet-50	SVM	88.1	90.0	92.8
		softmax	88.1	90.5	92.7

Table 15. TABLE XV: Comparison of iSQRT-COV and Improved B-CNN  [ 24 ] in terms of memory cost, training speed and classification accuracy (%). Following the experimental settings of Improved B-CNN  [ 24 ] , we employ VGG-VD16 as backbone model and set batch-size to 24, while performing fine-tuning on Birds-CUB200-2011 dataset. Here, we give memory taken by GCP layer for the whole mini-batch, percentage of memory overhead due to caching of intermediate tensors ( + Δ Δ +\Delta ) with respect to that of the whole networks (ALL), training speed (Hz) and classification accuracy. T indicates number of iterations.

Methods	FP	BP	Memory for GCP layer	$+\Delta$/ALL	Training speed (Hz)	Accuracy
Improved B-CNN [24]	NS iteration (T=5)	SVD+lyap	224M + 0M	0	17.13	85.8 [24]
	NS iteration (T=5)	iterative Lyap solver (T=5)		0	29.41	84.8
	NS iteration (T=10)	iterative Lyap solver (T=10)		0	28.30	85.2
	NS iteration (T=20)	iterative Lyap solver (T=20)		0	26.12	85.7
iSQRT-COV	NS iteration (T=5)	NS iteration (T=5)	224M + 96.18M	0.89%	31.25	87.2

Equations

$$\mathbf{\Sigma }=\mathbf{XJ}{\mathbf{X}}^T,\mathbf{J}=\frac{1}{M}(\mathbf{I}-\frac{1}{M}11^T),$$

$$\mathbf{\Sigma }=\mathbf{U}\mathbf{\Lambda }{\mathbf{U}}^T,$$

$$\mathbf{Z}\stackrel{△}{=}{\mathbf{\Sigma }}^{\alpha }=\mathbf{U}\mathrm{diag}(f({\lambda }_1),\dots ,f({\lambda }_d)){\mathbf{U}}^T.$$

$$f({\lambda }_i)={\lambda }_i^{\alpha }.$$

$$\arg \underset{\mathbf{\Sigma }}{\min }\log |\mathbf{\Sigma }|+\mathrm{tr}({\mathbf{\Sigma }}^{-1}\mathbf{S}),$$

$$\arg \underset{\hat{\mathbf{\Sigma }}}{\min }\log |\hat{\mathbf{\Sigma }}|+\mathrm{tr}({\hat{\mathbf{\Sigma }}}^{-1}\mathbf{S})+\gamma D_{\text{vN}}(\mathbf{I},\hat{\mathbf{\Sigma }}),$$

$${\mathbf{\Sigma }}^{\frac{1}{2}}=\arg \underset{\hat{\mathbf{\Sigma }}}{\min }\log |\hat{\mathbf{\Sigma }}|+\mathrm{tr}({\hat{\mathbf{\Sigma }}}^{-1}\mathbf{S})+D_{\mathrm{vN}}(\mathbf{I},\hat{\mathbf{\Sigma }}),$$

$$\displaystyle f(\lambda_{i})=\left\{\begin{matrix}{\lambda_{i}^{\alpha}}\Big{/}{\lambda_{1}^{\alpha}}&\text{for MPN-COV+M-$\ell_{2}$}\\ {\lambda_{i}^{\alpha}}\Big{/}(\sum_{k}\lambda_{k}^{2\alpha})^{\frac{1}{2}}&\text{for MPN-COV+M-Fro}\end{matrix}\right.$$

$$\displaystyle\mathrm{tr}\Big{(}\Big{(}\frac{\partial l}{\partial\mathbf{U}}\Big{)}^{T}\mathrm{d}\mathbf{U}+\Big{(}\frac{\partial l}{\partial\boldsymbol{\Lambda}}\Big{)}^{T}\mathrm{d}\boldsymbol{\Lambda}\Big{)}=\mathrm{tr}\Big{(}\Big{(}\frac{\partial l}{\partial\mathbf{Z}}\Big{)}^{T}\mathrm{d}\mathbf{Z}\Big{)},$$

$$\frac{\partial l}{\partial \mathbf{U}}$$

$$\frac{\partial l}{\partial \mathbf{\Lambda }}$$

$$\frac{\partial l}{\partial \mathbf{\Lambda }}=$$

$$\displaystyle-\mathrm{diag}\bigg{(}\dfrac{\alpha}{\lambda_{1}}\mathrm{tr}\Big{(}\mathbf{Z}\dfrac{\partial l}{\partial\mathbf{Z}}\Big{)},0,\ldots,0\bigg{)}$$

$$\frac{\partial l}{\partial \mathbf{\Lambda }}=$$

$$-$$

$$\displaystyle\dfrac{\partial l}{\partial\boldsymbol{\Sigma}}=\mathbf{U}\Big{(}\Big{(}\mathbf{K}^{T}\circ\Big{(}\mathbf{U}^{T}\dfrac{\partial l}{\partial\mathbf{U}}\Big{)}\Big{)}+\Big{(}\dfrac{\partial l}{\boldsymbol{\partial\Lambda}}\Big{)}_{\mathrm{diag}}\Big{)}\mathbf{U}^{T},$$

$$\displaystyle\dfrac{\partial l}{\partial\mathbf{X}}=\bigg{(}\dfrac{\partial l}{\partial\boldsymbol{\Sigma}}+\bigg{(}\dfrac{\partial l}{\partial\boldsymbol{\Sigma}}\bigg{)}^{T}\bigg{)}\mathbf{X}\mathbf{J}.$$

$$\displaystyle p(\mathbf{x})=\frac{1}{(2\pi)^{\frac{d}{2}}|\boldsymbol{\Sigma}|^{\frac{1}{2}}}\exp\big{(}-\frac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{T}\boldsymbol{\Sigma}^{-1}(\mathbf{x}-\boldsymbol{\mu})\big{)},$$

$$\varphi :\mathcal{N}(\mathbf{\mu },\mathbf{\Sigma })\stackrel{{\mathbf{\Sigma }}^{-1}=\mathbf{L}{\mathbf{L}}^T}{\to }{\mathbf{H}}_{\mathbf{\mu },\mathbf{J}}\stackrel{△}{=}\left[\begin{array}{cc}\mathbf{J}& \mathbf{\mu }\\ 0^T& 1\end{array}\right],$$

$$\psi :{\mathbf{H}}_{\mathbf{\mu },\mathbf{J}}\stackrel{{\mathbf{H}}_{\mathbf{\mu },\mathbf{J}}={\mathbf{S}}_{\mathbf{\mu },\mathbf{J}}{\mathbf{Q}}_{\mathbf{\mu },\mathbf{J}}}{\to }{\mathbf{S}}_{\mathbf{\mu },\mathbf{J}},$$

$$\mathcal{N}(\mathbf{\mu },\mathbf{\Sigma })\stackrel{\psi \circ \varphi }{\to }{\mathbf{S}}_{\mathbf{\mu },\mathbf{J}}={\left[\begin{array}{cc}\mathbf{\Sigma }+\mathbf{\mu }{\mathbf{\mu }}^T& \mathbf{\mu }\\ {\mathbf{\mu }}^T& 1\end{array}\right]}^{\frac{1}{2}}.$$

$$\mathbf{Y}=\frac{1}{M}\mathbf{AX}{\mathbf{X}}^T{\mathbf{A}}^T+\frac{2}{M}{\left(AX1{\mathbf{b}}^T\right)}_{\mathrm{sym}}+\left[\begin{array}{cc}\mathbf{O}& 0\\ 0^T& 1\end{array}\right]$$

$$\mathbf{Y}=\tilde{\mathbf{U}}\tilde{\mathbf{\Lambda }}{\tilde{\mathbf{U}}}^T,$$

$$\mathbf{Z}=\tilde{\mathbf{U}}{\tilde{\mathbf{\Lambda }}}^{\frac{1}{2}}{\tilde{\mathbf{U}}}^T,$$

$$\displaystyle\frac{\partial l}{\partial\mathbf{X}}=\frac{2}{M}\mathbf{A}^{T}\Big{(}\frac{\partial l}{\partial\mathbf{Y}}\Big{)}_{\mathrm{sym}}\Big{(}\mathbf{A}\mathbf{X}+\mathbf{b}\mathbf{1}^{T}\Big{)}.$$

$${\mathbf{Y}}_k$$

$${\mathbf{P}}_k$$

$${\mathbf{Y}}_k$$

$${\mathbf{P}}_k$$

$$\mathbf{A}=\frac{1}{\mathrm{tr}(\mathbf{\Sigma })}\mathbf{\Sigma }\text{or}\frac{1}{\Vert \mathbf{\Sigma }{\Vert }_F}\mathbf{\Sigma }.$$

$$\mathbf{Z}=\sqrt{\mathrm{tr}(\mathbf{\Sigma })}{\mathbf{Y}}_N\text{or}\mathbf{Z}=\sqrt{\Vert \mathbf{\Sigma }{\Vert }_F}{\mathbf{Y}}_N.$$

$$\frac{\partial l}{\partial \mathbf{\Sigma }}=$$

$$\displaystyle+\dfrac{1}{2\sqrt{\mathrm{tr}(\boldsymbol{\Sigma)}}}\mathrm{tr}\Big{(}\Big{(}\dfrac{\partial l}{\partial\mathbf{Z}}\Big{)}^{T}\mathbf{Y}_{N}\Big{)}\mathbf{I}.$$

Full text

1 Introduction

Deep learning methods, particularly deep convolutional neural networks (CNNs), have achieved great success in many computer vision tasks, especially in visual recognition [1, 2, 3, 4, 5]. As described in [6], deep learning automatically learns composition of multiple transformations that aims to produce abstract representations for predictions, and a good representation is one that can capture distribution of input data and is discriminative as an input of supervised prediction. However, existing deep CNNs often use global average pooling (GAP) to summarize the last convolution features as input of prediction, which captures only first-order statistics, limiting representation and generalization abilities of deep CNNs. To overcome this problem, one possible solution is to replace first-order GAP by some more powerful statistical modeling methods.

Compared with first-order pooling methods, covariance (second-order) pooling is able to capture richer statistics of features so that it can generate more informative representations [7, 8, 9]. In the classical, shallow architectures, Tuzel et al. [7] for the first time utilize covariance matrices for representing regular regions of images. Carreira et al. [8] present an O2P method, where second-order non-central moments are used to model free-form regions of images. Since covariance matrices are symmetric positive definite (SPD) matrices whose space forms a Riemannian manifold, geometric structure should be favorably considered for realizing the full potential. Therefore, affine invariant Riemannian metric (AIRM) [10] and Log-Euclidean Riemannian metric (LERM) [11] are respectively employed in [7] and [8] to measure the distances between covariance matrices. In particular, many researches  [12, 8, 9, 13] show that covariance representations with hand-crafted features (e.g., SIFT [14]) significantly outperform first-order counterparts.

In deep architectures, Ionescu et al. propose DeepO2P networks [15] in which second-order pooling is inserted after the last convolution layer of deep CNNs. They establish theory and practice of backpropagation, which enables end-to-end learning of CNNs involving non-linear structural matrices. A parallel work is bilinear CNN (B-CNN) [16] model, which performs outer product pooling of output features of the last convolution layers from two CNNs, producing non-central moments when the two CNNs are identical. Their main difference lies in normalization on SPD matrices obtained by covariance pooling: DeepO2P uses matrix logarithm while B-CNN performs element-wise power normalization followed by $\ell_{2}$-normalization.

Although above global covariance pooling (GCP) methods [15, 16] have been studied in deep architectures and report improvement, there exist two challenges that remain unresolved. The first one is robust covariance estimation. The dimensions of output features of the last convolution layers in deep CNNs usually are very high (e.g., 256 in AlexNet [1], 512 in VGG-VD [2] and 2048 in ResNet [4]), but the feature number is small ($<200$). Under such a scenario of high dimension and small sample size (HDSS), it is well-known that the sample covariance obtained by the classical maximum likelihood estimation (MLE) is not robust [19, 20]. Furthermore, recent studies [9, 17] have shown that robust covariance estimation leads to clear performance improvement in the case of HDSS. The second challenge is how to make favorable use of Riemannian geometry of covariance matrices. DeepO2P [15] exploits geometry of covariances using LERM [11], but our experiments (Section 4.2.2) show that LERM brings side effect for covariance matrices in the HDSS scenario. B-CNN [16] simply regards the space of covariances as a Euclidean space, discarding its geometric structure. Furthermore, neither DeepO2P [15] nor B-CNN [16] is concerned with robust covariance estimation.

Therefore, a question naturally arises: Can we overcome above challenges that GCP faces to further improve the representation and generalization abilities of deep CNNs? Inspired by a regularized maximum likelihood estimator (namely vN-MLE) of our previous work [17] and Power-Euclidean (Power-E) metric [18], this paper proposes a global Matrix Power Normalized COVariance (MPN-COV) Pooling to address the above two challenges. As shown in [17], vN-MLE performs robustly in the HDSS scenario, superior to most of the robust covariance estimators. The solution of vN-MLE corresponds to shrinkage of eigenvalues of covariance matrices [21, 22], and, interestingly, one special solution of vN-MLE is matrix square root (power of 1/2) of sample covariance. On the other hand, Power-E metric [18], which computes matrix power of SPD matrices for measuring distances, has been successfully used in medical imaging. Power-E metric has close connection with LERM [11] and so can approximately measure the Riemannian distance on the space of covariances. As such, our matrix power normalization in the context of deep CNNs amounts to robust covariance estimation while approximately yet effectively exploiting Riemannian geometry of covariances. The statistical and geometrical mechanisms underlying MPN-COV and deep CNNs with MPN-COV (namely MPN-COV-Net) are presented in Section 3.1 and Section 3.2, respectively.

Although MPN-COV can handle aforementioned challenges, it has several downsides. Firstly, it discards first-order statistics (mean vector) that is widely used in conventional deep CNNs. In order to integrate first-order information, we use a Gaussian model that accommodates both mean vector and covariance matrix for further performance improvement. We insert the Gaussian, identified as matrix square root of an SPD matrix, into deep CNNs based on Lie group theory [23]. The resulting network, called G2DeNet, is presented in Section 3.3. Secondly, forward propagation (FP) and backward propagation (BP) of both MPN-COV and G2DeNet need to compute matrix power that depends on eigen-decomposition (EIG) or singular value decomposition (SVD). However, implementation of EIG/SVD is limitedly supported on CUDA platform, even much slower than their CPU counterparts [15, 24]), resulting in computational bottleneck. Inspired by [24], we develop an iterative matrix square root normalized covariance pooling (iSQRT-COV)111The matrix square root normalization is a special case of MPN-COV where the power is 1/2 that usually performs best among all values of power (Refer to experiments in Sec. 4.2.1). based on Newton-Schulz iteration [25] for end-to-end learning, where we introduce key pre-normalization and post-compensation for the Newton-Schulz iteration, without which deeper CNNs (e.g., ResNet) fail to converge. The iSQRT-COV, as described in Section 3.4, is very suitable for parallel implementation on GPU, which can significantly speed up training of the networks. Thirdly, the output features of deep CNNs usually are of high dimension, resulting in much larger size of covariance representations that lead to high model complexity. To reduce model complexity without sacrificing performance, in Section 3.5, we propose to exploit progressive $1\times 1$ convolutions and group convolution [26] to compress covariance representations. All these efforts bring about a unified methodology, in which we not only improve significantly the existing GCP but also achieve fast training speed and affordable model complexity, making our final GCP networks very competitive and appealing.

The proposed GCP methods can be readily inserted into existing deep CNN architectures in an end-to-end manner. Fig. 1 illustrates deep CNNs with our proposed MPN-COV and improved solutions. Finally, we conduct experiments on a variety of visual recognition tasks, including large-scale object classification, scene categorization, fine-grained visual recognition, and texture classification. This paper summarizes and extends our preliminary works [17, 27, 28, 29], and the contributions are summarized as follows:

• –

We propose a global Matrix Power Normalized COVariance (MPN-COV) pooling for deep CNNs, which can address the challenges of robust covariance estimation and usage of Riemannian geometry of covariances, further improving the representation and generalization abilities of deep CNNs.

• –

We propose several solutions to overcome the downsides of MPN-COV. First, we propose a Gaussian embedding network to properly incorporate additional first-order information. Then, we implement forward and backward propagations of MPN-COV with $\alpha=1/2$ based on Newton-Schulz iteration for fast training of MPN-COV networks. Additionally, a compact strategy of progressive $1\times 1$ convolutions and group convolution is introduced to reduce size of covariance representations. These solutions as well as MPN-COV constitute our matrix power normalization methodology for improving deep CNNs.

• –

The proposed methods are implemented on different deep learning platforms, and the complete code will be released to open source repository. We write C++ code based on NVIDIA cuBLAS and Matlab using MatConvNet [30]. Meanwhile, we implement iSQRT-COV using PyTorch and TensorFlow packages.

• –

Extensive evaluations are conducted on various visual recognition tasks, including large-scale object classification on ImageNet [31] and scene categorization on Places365 [32], fine-grained visual recognition on Birds-CUB200-2011 [33], Aircrafts [34] and Cars [35], texture classification on DTD [36] and Indoor67 [37], and iNaturalist Challenge 2018 held in FGVC5 workshop in conjunction with CVPR 2018. The results with different CNN architectures (e.g., AlexNet [1], VGG-VD [2], ResNet [4] and DenseNet [5]) show the superiority of our GCP networks.

2 Related Work

This section reviews the related works that improve deep CNNs by integration of trainable sophisticated pooling or encoding methods, which are divided into four categories.

2.1 Deep CNNs with Global Second-order Pooling

Both DeepO2P [15] and B-CNN [16] insert a trainable second-order non-central moment into deep CNNs, where matrix logarithm normalization and element-wise power normalization followed by $\ell_{2}$-normalization are performed, respectively. Acharya et al. [38] explore a manifold network structure of covariance pooling for facial expression recognition, which performs, for input covariance matrix, bilinear mapping for dimension reduction, eigenvalue rectification and matrix logarithm consecutively. In improved B-CNN [24], Lin and Maji study the effect of different normalization methods on second-order statistics, finding out matrix square root normalization offers significant improvement over other normalization methods. In particular, they for the first time propose to use Newton-Schulz iteration [25] for efficient FP while computing accurate gradients by Layapnov equation for BP. While the idea of matrix square root normalization (a special case of MPN with power of 1/2) on GCP in improved B-CNN shares similarity with our MPN-COV [27], we provide theoretical explanations on why matrix power normalization works in deep architectures and experiment with large-scale ImageNet. Our iSQRT-COV [29], which aims to speed up our MPN-COV, is inspired by improved B-CNN but has clear differences in several respects. First, we develop a meta-layer consisting of pre-normalization, coupled matrix iteration and post-compensation, where pre-normalization (by trace or Frobenius norm) and post-compensation are distinctly important for convergence of deep CNNs (e.g., ResNet-50), and have been not studied in previous works including improved B-CNN. Second, both FP and BP of iSQRT-COV are performed based on Newton-Schulz iteration, which is more efficient than improved B-CNN. Finally, our method outperforms improved B-CNN on both large-scale and small-scale classification, while making first attempt to show GCP can benefit deeper CNNs (e.g., ResNet-50).

2.2 Deep CNNs with Global Approximate High-order Pooling

In general, high-order pooling methods result in large-size representations. Gao et al. [39] and Kong et al. [40] propose compact B-CNN and low-rank B-CNN models to reduce sizes of covariance (second-order) representations, respectively. These methods replace exact covariances by small-size approximation ones, while achieving comparable performance. Based on compact B-CNN model [39], Dai et al. [41] fuse additional first-order (mean) information by simply concatenating them. Kernel pooling [42] extends the similar idea with the compact B-CNN to approximate higher-order (number of order $>2$) pooling. Sharing similar philosophy with [42], Cai et al. [43] obtain small-size higher-order representations based on polynomial kernel approximation and rank-$1$ tensor decomposition [44], namely HIHCA. Both kernel pooling [42] and HIHCA [43] improve compact B-CNN by exploiting higher-order information. In [45], Yu and Salzmann propose a statistically-motivated second-order (SMSO) pooling, which successively consists of parametric vectorization, an element-wise square-root normalization, and a trainable affine transformation, each of which yields a well-defined distribution of data. SMSO can produce compact second-order representation while achieving state-of-the-art results. However, all above methods consider neither robust estimation nor geometry of manifold, limiting the ability of higher-order statistical modeling.

2.3 Deep CNNs with Local Second-order Statistics

In contrary to above works exploring global second-order pooling, some researchers try to incorporate local second-order statistics into deep CNNs. Among them, Factorized Bilinear (FB) method [46] introduces an additional parametric quadratic term into linear transformation of convolution or fully-connected (FC) layers. FB can incorporate more complex non-linearity structures into deep CNNs by considering second-order interaction between information flow. Second-Order Response Transform (SORT) [47] proposes to fuse outputs of two-branch block using a second-order term (i.e., element-wise product and a sum operation), in order to increase the nonlinearity of deep CNNs as well. Obviously, these methods also discard robust estimation and geometry of second-order statistics so that they cannot make full use of capacity of second-order pooling.

2.4 Deep CNNs with Trainable BoVW Methods

In the past decades, Bag-of-Visual-Words (BoVW) model is one of the most widely used orderless pooling methods for visual classification. Recently, some works insert high-performance BoVW methods [48, 12] as trainable structural layers into deep CNNs. Thereinto, NetVLAD [49] implements the modified vector of locally aggregated descriptors (VLAD) [48] in an end-to-end manner. FisherNet [50] accomplishes the trainable layer of simplified Fisher vector (FV) [12]. Unlike FisherNet, Li et al. [51] propose a MFAFVNet for scene categorization, which performs a deep embedded implementation of mixture of factor analyzers Fisher vector (MFA-FV) method [52]. Different from these methods, this paper proposes to integrate a global matrix power normalized covariance pooling into deep CNNs.

3 The Proposed Method

In this section, we first describe our MPN-COV and the underlying mechanisms, and then instantiate our MPN-COV-Net. Subsequently, we present a global Gaussian embedding network to fuse first-order information, develop an iterative matrix square root normalization method to speed up training of networks, and introduce a compact strategy to reduce size of covariance representations.

3.1 Global Matrix Power Normalized COVariance (MPN-COV) Pooling

3.1.1 Computation of MPN-COV

Let $\mathcal{X}\in\mathbb{R}^{w\times h\times d}$ be the outputs of the convolution layer right before computation of MPN-COV, where $w$, $h$ and $d$ indicate spatial width, height and the number of channels, respectively. Then, the feature tensor $\mathcal{X}$ is reshaped to a feature matrix $\mathbf{X}\in\mathbb{R}^{d\times M}$ where $M=w\times h$. Given the feature matrix $\mathbf{X}$ consisting of $M$ samples with $d$-dimension, its sample covariance is computed as

[TABLE]

where $\mathbf{I}$ indicates a $M\times M$ identity matrix, $\mathbf{1}$ is a $M$-dimension vector with all elements being one, and $T$ denotes the matrix transpose.

Sample covariance $\boldsymbol{\Sigma}$ is a symmetric positive definite or semidefinite matrix, which can be factorized by EIG/SVD:

[TABLE]

where $\boldsymbol{\Lambda}=\mathrm{diag}(\lambda_{1},\ldots,\lambda_{d})$ is a diagonal matrix and $\lambda_{i},i=1,\ldots,d$ are eigenvalues arranged in non-increasing order; $\mathbf{U}=[\mathbf{u}_{1},\ldots,\mathbf{u}_{d}]$ is an orthogonal matrix whose column $\mathbf{u}_{i}$ is the eigenvector corresponding to $\lambda_{i}$. Through EIG or SVD, we can compute matrix power as follows:

[TABLE]

Here $\alpha>0$ is a scalar and $f(\lambda_{i})$ is power of the eigenvalues

[TABLE]

In this paper, the operation in Eq. (3) is called MPN-COV, which is inserted after the last convolution layer of deep CNNs to collect second-order statistics of the output features (i.e., $\mathbf{X}$) as global representations. Next, we describe the mechanisms underlying MPN-COV in terms of robust covariance estimation and usage of geometry.

3.1.2 Robust Covariance Estimation

Assuming $\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{M}]$ are sampled from a Gaussian distribution, the covariance $\boldsymbol{\Sigma}$ of $\mathbf{X}$ can be estimated by optimizing the following objective function based on MLE:

[TABLE]

where $\mathbf{S}=\mathbf{X}\mathbf{J}\mathbf{X}^{T}$ is sample covariance, $|\cdot|$ indicates matrix determinant and $\mathrm{tr}(\cdot)$ means trace of matrix. The solution to the MLE (5) is the sample covariance, i.e., $\boldsymbol{\Sigma}=\mathbf{S}$. However, it is well known that MLE is not robust when data is of HDSS [19, 20]. This scenario is just what our covariance pooling faces: in most of the existing deep CNNs [2, 3, 4, 5], the output features (i.e., $\mathbf{X}$) of last convolution layer have less than 200 samples of dimension larger than 256, so the sample covariances are always ill-conditioned, rendering robust estimation critical.

For robust estimation of covariances in the case of HDSS, the general principle is shrinkage of eigenvalues of the sample covariances for counteracting the ill-conditioning of covariance matrices [21, 22]. Besides, some researchers propose various regularized MLE methods for robust covariance estimation (see [53] and references therein). Notably, our MPN-COV closely conforms to the shrinkage principle [21, 22], i.e., shrinking the largest sample eigenvalues and stretching the smallest ones. Moreover, MPN-COV can be deemed as a regularized MLE, namely vN-MLE  [17]:

[TABLE]

where $\gamma>0$ is a regularizing constant, and $D_{\text{vN}}(\mathbf{A},\mathbf{B})=\mathrm{tr}(\mathbf{A}(\log(\mathbf{A})-\log(\mathbf{B}))-\mathbf{A}+\mathbf{B})$ is matrix von-Neumann divergence. Compared with the classical MLE (5) which only includes the first two terms, the objective function of vN-MLE (6) introduces the third term, constraining the estimated covariance $\widehat{\boldsymbol{\Sigma}}$ be similar to the identity matrix $\mathbf{I}$. In [17], it has been shown that the vN-MLE outperforms other shrinkage methods [21, 22, 54] and regularized MLE method [53]. Briefly, we have

Proposition 1

MPN-COV with $\alpha=\frac{1}{2}$ is the unique solution to the vN-MLE in which $\gamma=1$, i.e.,

[TABLE]

where $\boldsymbol{\Sigma}=\mathbf{S}$.

Proposition 1 shows our MPN-COV with $\alpha=1/2$ performs robust estimation of covariance, and experiments in Section 4.2.1 show that $\alpha=1/2$ performs best. Details on vN-MLE can be referred to [17].

3.1.3 Approximate Usage of Geometry

Since the space of $d\times d$ covariance matrices (denoted by $\mathbb{S}^{+}_{d}$) forms a Riemannian manifold, geometry should be considered when distances between covariances are measured. There are mainly two kinds of Riemannian metrics, i.e., AIRM [10] and LREM [11]. The AIRM is affine-invariant, but it is computationally inefficient and coupled, not scalable to large-scale scenarios. The LERM is a similarity-invariant decoupled metric and efficient to compute, so that it is scalable to large-scale problems. Our MPN-COV can be regarded as matching covariance matrices with the Power-Euclidean (Pow-E) metric [18], which has close connection with the LERM, as presented in the following proposition:

Proposition 2

For any two covariance matrices $\boldsymbol{\Sigma}_{1}$ and $\boldsymbol{\Sigma}_{2}$, the limit of the Pow-E metric d_{\alpha}(\boldsymbol{\Sigma}_{1},\boldsymbol{\Sigma}_{2})=\frac{1}{\alpha}\big{\|}\boldsymbol{\Sigma}^{\alpha}_{1}-\boldsymbol{\Sigma}_{2}^{\alpha}\big{\|}_{F} as $\alpha>0$ approaches 0 equals the LERM, i.e., \lim\limits_{\alpha\rightarrow 0}d_{\alpha}(\boldsymbol{\Sigma}_{1},\boldsymbol{\Sigma}_{2})=\big{\|}\log(\boldsymbol{\Sigma}_{1})-\log(\boldsymbol{\Sigma}_{2})\big{\|}_{F}.

This conclusion was first mentioned in [18] but without proof. Here we briefly prove this proposition. Note that d_{\alpha}(\boldsymbol{\Sigma}_{1},\boldsymbol{\Sigma}_{2})=\big{\|}\frac{1}{\alpha}(\boldsymbol{\Sigma}^{\alpha}_{1}-\mathbf{I})-\frac{1}{\alpha}(\boldsymbol{\Sigma}^{\alpha}_{2}-\mathbf{I})\big{\|}_{F}. For any covariance $\boldsymbol{\Sigma}$ we have $\frac{1}{\alpha}(\boldsymbol{\Sigma}^{\alpha}-\mathbf{I})=\mathbf{U}\mathrm{diag}(\frac{\lambda_{1}^{\alpha}-1}{\alpha},\ldots,\frac{\lambda_{n}^{\alpha}-1}{\alpha})\mathbf{U}^{T}$ based on its EIG. The identity about the limit in Proposition 2 follows immediately by recalling $\lim_{\alpha\rightarrow 0}\frac{\lambda^{\alpha}-1}{\alpha}=\log(\lambda)$. Hence, the proposed MPN-COV can be viewed as approximately exploiting the Riemannian geometry of $\mathbb{S}^{+}_{d}$. It seems that the LERM is better than the Pow-E metric, since the former computes the true geodesic distance but the latter only measures it approximately. We argue that this is not the case for the scenario of deep CNNs from the perspectives of both numerical stability and distribution of eigenvalues. Detailed discussion is given in Section 4.2.2.

3.2 Global MPN-COV Pooling Network

We first instantiate a global matrix power normalized covariance pooling neural network (MPN-COV-Net) by inserting our MPN-COV after the last convolution layer of deep CNNs, in place of the common GAP. The diagram of our MPN-COV block is illustrated in Fig. 2. For forward propagation, we first compute the covariance pooling of output features $\mathbf{X}$ of the last convolution layer using Eq. (1), then we perform matrix power normalization (3). Inspired by the element-wise power normalization technique [12], we can further perform, after MPN-COV, normalization by matrix $\ell_{2}-$norm (M-$\ell_{2}$) or by matrix Frobenius norm (M-Fro). The matrix $\ell_{2}-$norm (also known as the spectral norm) of a matrix $\boldsymbol{\Sigma}$, denoted by $\|\boldsymbol{\Sigma}\|_{2}$, is defined as the largest singular value of $\boldsymbol{\Sigma}$, which equals the largest eigenvalue if $\boldsymbol{\Sigma}$ is a covariance matrix. The matrix Frobenius norm of $\boldsymbol{\Sigma}$ can be defined in various ways such as $\|\boldsymbol{\Sigma}\|_{F}=(\mathrm{tr}(\boldsymbol{\Sigma}^{T}\boldsymbol{\Sigma}))^{\frac{1}{2}}=(\sum_{i}\lambda_{i}^{2})^{\frac{1}{2}}$, where $\lambda_{i}$ are eigenvalues of $\boldsymbol{\Sigma}$. Then, we have

[TABLE]

For backpropagation of MPN-COV block, we need to compute the partial derivative of loss function $l$ with respect to the input $\mathbf{X}$ based on the methodology of matrix backpropagation [15, 55]. First of all, given $\frac{\partial l}{\partial\mathbf{Z}}$ propagated from top layer, we compute the derivatives $\frac{\partial l}{\partial\mathbf{U}}$ and $\frac{\partial l}{\partial\boldsymbol{\Lambda}}$ based on the following chain rule:

[TABLE]

where $\mathrm{d}\mathbf{Z}$ denotes variation of matrix $\mathbf{Z}$. According to Eq.  (3), we have $\mathrm{d}\mathbf{Z}=\mathrm{d}\mathbf{U}\mathbf{F}\mathbf{U}^{T}+\mathbf{U}\mathrm{d}\mathbf{F}\mathbf{U}^{T}+\mathbf{U}\mathbf{F}\mathrm{d}\mathbf{U}^{T}$, where $\mathbf{F}=\mathrm{diag}(f(\lambda_{1}),\ldots,f(\lambda_{d}))$ and \mathrm{d}\mathbf{F}=\mathrm{diag}\big{(}\alpha\lambda_{1}^{\alpha-1},\ldots,\alpha\lambda_{d}^{\alpha-1}\big{)}\mathrm{d}\boldsymbol{\Lambda}. After some arrangements, we obtain

[TABLE]

where $(\cdot)_{\mathrm{diag}}$ denotes the matrix diagonalization. For MPN-COV+M-$\ell_{2}$ and MPN-COV+M-Fro, $\frac{\partial l}{\partial\boldsymbol{\Lambda}}$ takes respectively the following forms:

[TABLE]

and

[TABLE]

Next, given $\frac{\partial l}{\partial\mathbf{U}}$ and $\frac{\partial l}{\partial\boldsymbol{\Lambda}}$, we need to compute $\frac{\partial l}{\partial\boldsymbol{\Sigma}}$ associated with Eq. (2), whose corresponding chain rule is $\mathrm{tr}((\frac{\partial l}{\partial\boldsymbol{\Sigma}})^{T}\mathrm{d}\boldsymbol{\Sigma})=\mathrm{tr}((\frac{\partial l}{\partial\mathbf{U}})^{T}\mathrm{d}\mathbf{U}+(\frac{\partial l}{\partial\boldsymbol{\Lambda}})^{T}\mathrm{d}\boldsymbol{\Lambda})$. Note that $\mathbf{U}$ is an orthogonal matrix. After some arrangements, we have

[TABLE]

where $\circ$ denotes matrix Hadamard product. The matrix $\mathbf{K}=\{K_{ij}\}$ where $K_{ij}=1/(\lambda_{i}-\lambda_{j})$ if $i\neq j$ and $K_{ij}=0$ otherwise. We refer readers to  [55, Proposition 2] for in-depth derivation of Eq. (13).

Finally, given $\frac{\partial l}{\partial\boldsymbol{\Sigma}}$, we derive the gradient of the loss function with respect to the input feature $\mathbf{X}$, which has

[TABLE]

As described above, by using the formulas of Eq. (3) and Eq. (14), our MPN-COV-Net can be end-to-end trained.

3.3 Global Gaussian Embedding Network

Our MPN-COV can address the two challenges GCP faces, clearly outperforming other competing methods. Nevertheless, it neglects the first-order representation (i.e., mean vector) that is widely used in conventional deep CNNs. Previous researches [12, 9, 17] have shown that combination of first- and second-order statistics is often better than either one single statistics. As is commonly done, we use Gaussians for modeling distributions of features. However, as the space of Gaussians is a manifold, it is challenging to insert a Gaussian into deep CNNs. As such, we propose to use Gaussian embedding method proposed in [23]. This embedding method equips the space of Gaussians with a Lie group structure, respecting both geometrical and algebraic structures, and meanwhile it identifies a Gaussian as matrix square root of an SPD matrix, suitable for backpropagation while achieving better performance. In Appendix II, we compare with other embedding methods [56, 57, 58, 59, 23]. In the following, we first describe Gaussian embedding method and then present FP and BP of the resulting network (i.e., G2DeNet).

Gaussian Embedding Let $\mathbf{X}=[\mathbf{x}_{1},\ldots,\mathbf{x}_{M}]\in\mathbb{R}^{d\times M}$ be a set of deep convolution features, we use Gaussian to model feature distribution as

[TABLE]

where $\boldsymbol{\mu}=\frac{1}{M}\sum_{i=1}^{M}\mathbf{x}_{i}$ and $\boldsymbol{\Sigma}=\frac{1}{M}\sum_{i=1}^{M}(\mathbf{x}_{i}-\boldsymbol{\mu})(\mathbf{x}_{i}-\boldsymbol{\mu})^{T}$ are mean vector and sample covariance matrix, respectively. Let $\boldsymbol{\Sigma}^{-1}=\mathbf{L}\mathbf{L}^{T}$ be the Cholesky decomposition of the inverse of $\boldsymbol{\Sigma}$, where $\mathbf{L}$ is a lower triangular matrix of order $d$ with positive diagonals. The Gaussian $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ can be uniquely mapped to a positive definite upper triangular matrices of order $d+1$ through the mapping $\phi$:

[TABLE]

where $\mathbf{J}=\mathbf{L}^{-T}$ and $\mathbf{H}_{\boldsymbol{\mu},\mathbf{J}}\in UT^{+}(d+1)$. $UT^{+}(d+1)$ indicates the set of all positive definite upper triangular matrices of order $d+1$. We note that the embedding (15) is not suitable for backpropagation due to Cholesky decomposition and matrix inverse. The matrix $\mathbf{H}_{\boldsymbol{\mu},\mathbf{J}}$ can be further mapped to a unique SPD matrix through a mapping $\psi$ based on its matrix polar decomposition $\mathbf{H}_{\boldsymbol{\mu},\mathbf{J}}=\mathbf{S}_{\boldsymbol{\mu},\mathbf{J}}\mathbf{Q}_{\boldsymbol{\mu},\mathbf{J}}$, i.e.,

[TABLE]

where $\mathbf{S}_{\boldsymbol{\mu},\mathbf{J}}$ is an SPD matrix and $\mathbf{Q}_{\boldsymbol{\mu},\mathbf{J}}$ is the closest orthogonal matrix to $\mathbf{H}_{\boldsymbol{\mu},\mathbf{J}}$. Through the two consecutive mappings $\phi$ and $\psi$, we identify a Gaussian as matrix square root of an SPD matrix:

[TABLE]

We suggest readers refer to [23] for theoretical details.

FP and BP of G2DeNet According to the embedding form in Eq. (17), we design global Gaussian embedding block which consists of two layers, i.e., matrix partition layer and square root SPD matrix layer, as shown in Fig. 3. Below we first describe the matrix partition layer. We let $\mathbf{Y}=\mathbf{S}_{\boldsymbol{\mu},\mathbf{J}}^{2}$, in which the mean vector $\boldsymbol{\mu}$ and covariance matrix $\boldsymbol{\Sigma}$ are obviously entangled. The purpose of this layer is to decouple $\mathbf{Y}$ and explicitly write it as the function of input feature matrix $\mathbf{X}$ to facilitate computation of the derivatives. We note that there exists the identity $\boldsymbol{\Sigma}=\frac{1}{M}\mathbf{X}\mathbf{X}^{T}-\boldsymbol{\mu}\boldsymbol{\mu}^{T}$. After some elementary manipulations, we have

[TABLE]

where $\mathbf{A}=\begin{bmatrix}\mathbf{I}&\mathbf{0}\end{bmatrix}^{T}$ consists of a $d\times d$ identity matrix $\mathbf{I}$ and a $d-$dimensional zero vector $\mathbf{0}$, $\mathbf{b}$ is a $(d+1)-$dimensional vector with all elements being zero except the last one which is equal to one, $\mathbf{1}$ is a $M-$dimensional vector with all elements being one, and $\mathbf{O}$ is a $d\times d$ zero matrix. The notation $(\mathbf{F})_{\mathrm{sym}}=\frac{1}{2}\left(\mathbf{F}+\mathbf{F}^{T}\right)$ denotes matrix symmetrization.

The purpose of Square Root SPD Matrix Layer is to compute the matrix square root of SPD matrix $\mathbf{Y}$, i.e., $\mathbf{Z}=\mathbf{Y}^{\frac{1}{2}}$. Similar to matrix power as described previously, the matrix $\mathbf{Z}$ can be computed via EIG/SVD. It is straightforward to know that $\mathbf{Y}$ is an SPD matrix, and so it can be factorized as

[TABLE]

where $\widetilde{\boldsymbol{\Lambda}}=\mathrm{diag}(\widetilde{\lambda}_{1},\cdots,\widetilde{\lambda}_{d+1})$ and $\widetilde{\mathbf{U}}=[\widetilde{\mathbf{u}}_{1}\;\cdots\;\widetilde{\mathbf{u}}_{d+1}]$ are the eigenvalues and eigenvectors of $\mathbf{Y}$, respectively. The matrix square root of $\mathbf{Y}$ can be computed as:

[TABLE]

where $\widetilde{\boldsymbol{\Lambda}}^{\frac{1}{2}}=\mathrm{diag}(\widetilde{\lambda}_{1}^{\frac{1}{2}},\cdots,\widetilde{\lambda}_{d+1}^{\frac{1}{2}})$ is computed as element-wise square root of the eigenvalues. Combining matrix partition layer (18) with square root SPD matrix layer (20), we can accomplish FP of global Gaussian embedding block.

Given $\frac{\partial l}{\partial\mathbf{Z}}$ propagated from top layer, we perform BP by computing $\frac{\partial l}{\partial\mathbf{Y}}$ and $\frac{\partial l}{\partial\mathbf{X}}$ associated with square root SPD matrix layer (20) and matrix partition layer (18), respectively. Note that derivation of $\frac{\partial l}{\partial\mathbf{Y}}$ shares similar philosophy with that of Eq. (13), and we omit it for simplicity. Given $\frac{\partial l}{\partial\mathbf{Y}}$, we can compute $\frac{\partial l}{\partial\mathbf{X}}$ based on Eq. (18), which takes the following form:

[TABLE]

3.4 Iterative Matrix Square Root Normalized Covariance Pooling Network

In our MPN-COV-Net and G2DeNet, computation of matrix square root heavily depends on EIG or SVD. However, fast implementation of EIG or SVD on GPU is still an open problem. Table I presents the running time of EIG and SVD of a $256\times 256$ covariance matrix. Matlab (M) built-in CPU functions deliver over 10x and 5x speedups over the CUDA counterparts and built-in GPU functions, respectively. As such, both our MPN-COV-Net and G2DeNet opt for EIG or SVD on CPU for computing matrix square root, greatly restricting the training speed of the networks on GPUs.

To overcome this limitation, inspired by [24], we develop a fast training method, which is called iterative matrix square root normalization of covariance (iSQRT-COV) pooling, making use of Newton-Schulz iteration [25] in both FP and BP. At the core of iSQRT-COV is a meta-layer with loop-embedded directed graph structure, which consists of three consecutive structured layers, performing pre-normalization, coupled matrix iteration and post-compensation, respectively. The iSQRT-COV block is illustrated in Fig. 4, and the details are described as follows.

Newton-Schulz Iteration Higham [25] studied a class of methods for iteratively computing matrix square root. These methods, termed as Newton-Padé iterations, are developed based on the connection between matrix sign function and matrix square root, together with rational Padé approximation. Specifically, for computing the matrix square root of $\mathbf{A}$, given $\mathbf{Y}_{0}=\mathbf{A}$ and $\mathbf{P}_{0}=\mathbf{I}$, for $k=1,\cdots,N$, the coupled iteration takes the following form [25, Chap. 6.7]:

[TABLE]

where $p_{lm}$ and $q_{lm}$ are polynomials, and $l$ and $m$ are non-negative integers. According to Eq. (3.4), $\mathbf{Y}_{k}$ and $\mathbf{P}_{k}$ quadratically converge to $\mathbf{A}^{1/2}$ and $\mathbf{A}^{-1/2}$, respectively. However, it converges only locally, i.e., it converges if $\|\mathbf{I}-\mathbf{A}\|<1$ where $\|\cdot\|$ denotes any induced (or consistent) matrix norm. The family of coupled iteration is stable in that small errors in the previous iteration will not be amplified. The case of $l=0,m=1$ called Newton-Schulz iteration fits for our purpose as no GPU unfriendly matrix inverse is involved:

[TABLE]

where $k=1,\ldots,N$. Clearly Eq. (3.4) involves only matrix product, suitable for parallel implementation on GPU. Compared to accurate square root computed by EIG, one can obtain approximate solution with a small number of iterations $N$, which is determined by cross-validation.

Pre-normalization and Post-compensation As Newton-Schulz iteration only converges locally, we pre-normalize covariance $\boldsymbol{\Sigma}$ by its trace or Frobenius norm, i.e.,

[TABLE]

Let $\lambda_{i}$ be eigenvalues of $\boldsymbol{\Sigma}$, arranged in nondecreasing order. As $\mathrm{tr}(\boldsymbol{\Sigma})=\sum_{i}\lambda_{i}$ and $\|\boldsymbol{\Sigma}\|_{F}=\sqrt{\sum_{i}\lambda_{i}^{2}}$, it is easy to see that $\|\mathbf{I}-\mathbf{A}\|_{2}$, which equals to the largest singular value of $\mathbf{I}-\mathbf{A}$, is $1-\frac{\lambda_{1}}{\sum_{i}\lambda_{i}}$ and $1-\frac{\lambda_{1}}{\sqrt{\sum_{i}\lambda_{i}^{2}}}$ for the case of trace and Frobenius norm, respectively, both less than 1. Hence, the convergence condition is satisfied.

The above pre-normalization for covariance pooling nontrivially changes the data magnitudes, which produces side effect on network such that the prevalent ResNet [4] fails to converge. To counteract this adverse influence, after the Newton-Schulz iteration, we accordingly perform post-compensation as follows:

[TABLE]

Alternatively, one may consider Batch Normalization (BN) [60], which does work, successfully helping ResNet convergence. Compared to BN, our post-compensation achieves about 1% lower top-1 error with ResNet-50 architecture on ImageNet (see Section 4.2.1 for details).

Note that derivation of BP of iSQRT-COV is not straightforward. Despite autograd toolkits provided by some deep learning frameworks can accomplish this task automatically, the involved BP is still in a black box. Meanwhile, autograd sometimes brings uncertainty, e.g., autograd of PyTorch 0.3.0 or below cannot compute gradients of iSQRT-COV correctly. To make iSQRT-COV be self-contained and enable its implementation to be accessible when autograd toolkit is unavailable (e.g., the well-known Caffe and early versions of MatConvNet), we show the gradient of the loss function $l$ with respect to $\boldsymbol{\Sigma}$ and give the derivation details in the Appendix I for conciseness. Specifically, given $\frac{\partial l}{\partial\mathbf{Z}}$, we can compute the gradient for pre-normalization by trace as

[TABLE]

If pre-normalization by Frobenius norm is adopted, we have

[TABLE]

where $\frac{\partial l}{\partial\mathbf{A}}$ can be found in Appendix I. Finally, given $\frac{\partial l}{\partial\boldsymbol{\Sigma}}$, the gradient of $l$ with respect to $\mathbf{X}$ is given in Eq. (14).

3.5 Compact Covariance Representations

Given an input feature matrix $\mathbf{X}\in\mathbb{R}^{d\times M}$ consisting of $M$ samples with $d$-dimension, size of covariance representation is $d\times(d+1)/2$ after matrix vectorization by considering the symmetry. Taking $d=512$ or $d=1024$ as an example, the size of covariance representation is $131\text{,}328$ or $524\text{,}800$. Such large-size representation leads to high model complexity. To deal with this problem, we present a compact strategy for compressing covariance representations, as shown in Fig. 5. The compact strategy consists of progressive $1\times 1$ convolutions and group convolution. The formers are used to reduce feature dimension right before computation of covariance matrices, while the latter is employed to reduce size of covariance representations after GCP block.

Since the dimension of input features decides size of covariance representation, we first introduce a strategy of progressive $1\times 1$ convolutions to reduce the dimension of input features. The $1\times 1$ convolution is first introduced in [61], and has been widely used in deep CNN architectures [3, 4, 5] for parameter reduction. However, an abrupt dimensionality reduction (e.g., directly from $2080$ to $128$) may hurt the performance of covariance representations. Therefore, we propose to progressively reduce the dimension of input features from $d$ to $\hat{d}$ ($\hat{d}\ll d$) using a set of consecutive $1\times 1$ convolutions, i.e., $d=d_{0}\rightarrow d_{1}\rightarrow\dots\rightarrow d_{K}=\hat{d}$ with $d_{k-1}>d_{k}$. In this way, dimension of covariance representation can be effectively reduced to $\hat{d}\times(\hat{d}+1)/2$ (e.g., $8256$ for $\hat{d}=128$). Compared with dimensionality reduction using only one $1\times 1$ convolution, our progressive convolutions can achieve better performance.

Once dimension of input is determined, we can further reduce the size of covariance representation by learning a low-dimensional linear transform matrix $\mathbf{W}$ on covariance representation $\mathbf{z}$, which corresponds to inserting a FC layer ($\mathbf{\hat{z}}=\mathbf{W}\mathbf{z}+\mathbf{b}$) after covariance pooling. The number of parameters of $\mathbf{W}$ is $d_{\mathbf{z}}\times d_{\mathbf{\hat{z}}}$, where $d_{\mathbf{z}}$ and $d_{\mathbf{\hat{z}}}$ are dimensions of $\mathbf{z}$ and $\mathbf{\hat{z}}$, respectively. If $d_{\mathbf{z}}$ is very high, transform matrix $\mathbf{W}$ will involve of many parameters, increasing computation and memory costs. Therefore, we introduce group convolution [26] to alleviate this problem. Different from standard convolution, group convolution divides $\mathbf{z}$ into $G$ groups and performs convolution operation in each sub-group. As such, number of parameters of $\mathbf{W}$ can be decreased to $(d_{\mathbf{z}}\times d_{\mathbf{\hat{z}}})/G$, very suitable for large-size covariances.

4 Experiments

To evaluate the representation and generalization abilities of the proposed methods, we conduct experiments on both large-scale object classification and scene categorization, small-scale fine-grained visual recognition and texture classification as well as large-scale species classification competition. We first describe implementation details of our methods. Then we make ablation study of the proposed methods on large-scale ImageNet dataset, and compare with state-of-the-art methods on ImageNet [31] and Places365 [32]. Additionally, we verify the generalization ability of the proposed methods through transferring them to three fine-grained (i.e., Birds-CUB200-2011 [33], Aircrafts [34] and Cars [35]) and two texture image datasets (i.e., Indoor67 [37] and DTD [36]). Finally, we show the results on iNaturalist Challenge 2018 held in CVPR 2018 workshop on FGVC5.

4.1 Implementation Details

To implement MPN-COV block, we use the EIG algorithm on CPU in single-precision floating-point format, as the EIG algorithm on GPU is much slower, despite the overhead due to data transfer between GPU memory and CPU memory. Since MPN-COV allows non-negative eigenvalues, we truncate to zeros the eigenvalues smaller than $\mathrm{eps}(\lambda_{1})$, which is a matlab function denoting the positive distance from the largest eigenvalue $\lambda_{1}$ to its next larger floating-point number. For achieving the global Gaussian embedding block, as suggested in [15], we use SVD to compute square root SPD matrix because SVD is numerically more stable, and SVD algorithm is implemented on CPU as well. Meanwhile, we add a small positive number 1e-3 to the diagonal entries of Gaussian embedded matrices for numerical stability. The implementation of iSQRT-COV block is encapsulated in three computational modules, which accomplish forward&backward computation of pre-normalization layer, Newton-Schulz iteration layer and post-compensation layer, respectively. We implement above methods based on the MatConvNet package [30]; iSQRT-COV-Net is also implemented using Pytorch and TensorFlow.

The proposed GCP blocks are inserted after the last convolution layer (with ReLU) of deep CNNs. Unless otherwise stated, we discard the last downsampling in the networks so that we have larger number of features, and add one $1\times 1$ convolution with $d=256$ channels before our covariance pooling for all CNN architectures but AlexNet. As such, we have $w\times h\times 256$ feature maps right before covariance pooling, where $w$ and $h$ are width and height of feature maps, obtaining a $256\times 256$ covariance matrix (or 32,896-$d$ covariance representation after matrix vectorization). We train our networks using mini-batch stochastic gradient descent algorithm with a momentum of 0.9 and a weight decay of $5\times 10^{-4}$. All programs run on two PCs each of which is equipped with a 4-core Intel i7-4790k@4.0GHz CPU, 32G RAM, 512GB Samsung PRO SSD and two Titan Xp GPUs. The remaining hyper-parameters of training networks will be described where appropriate in following subsections.

4.2 Large-scale Object Classification on ImageNet

4.2.1 Ablation Study

We first analyze the proposed methods on large-scale ImageNet dataset, which contains 1.28M training images, 50K validation images and 100K testing images collected from 1K classes. On this dataset, we follow [2] for data augmentation, and adopt the commonly used 1-crop or 10-crop prediction for performance evaluation. Following the common settings [1, 2, 4], we report the results on the validation set. Specifically, we study the effect of parameter $\alpha$ on MPN-COV and compare various normalization methods on GCP, as well as assess impacts of number $N$ of Newton-Schulz iterations and post-compensation on iSQRT-COV, speedup ratio of training iSQRT-COV-Net over MPN-COV-Net, and effect of compact covariance representations.

MPN-COV: Effect of parameter $\alpha$ Our MPN-COV has a key parameter $\alpha$, and we evaluate it under AlexNet architecture. Our MPN-COV-Net is trained up to 20 epochs, where the learning rates follow exponential decay, changing from $10^{-1.2}$ to $10^{-5}$ with a batch size of 128. Fig. 6 shows top-1 error vs. $\alpha$ using single-crop prediction. We first note that the Plain-COV ($\alpha=1$, no normalization) produces an error rate of 40.41%, about $1.1\%$ less than the original AlexNet. When $\alpha<1$, the normalization function shrinks eigenvalues larger than 1.0 and stretches those less than 1.0. As $\alpha$ (less than 1.0) decreases, the error rate continuously gets smaller until the smallest value at around $\alpha=\frac{1}{2}$. With further decline of $\alpha$, however, we observe the error rate grows consistently and soon is larger than that of the Plain-COV. Note that over the interval $[0.4,0.9]$ the performance of MPN-COV varies insignificantly. When $\alpha>1$, the effect of normalization is contrary, i.e., eigenvalues less than 1.0 are shrunk while those larger than 1.0 are stretched, which is not beneficial for covariance representations as indicated by the consistent growth of the error rates. We note that similar trend of effect of matrix power on classification performance was also reported in [24]. As such, we set $\alpha=\frac{1}{2}$ throughout the following experiments.

MPN-COV: Comparison of various normalization methods We compare four kinds of normalization methods, i.e., MPN, M-Fro, M-$\ell_{2}$ and element-wise power normalization followed by $\ell_{2}$-normalization (E-PN for short) [12, 16]. Table II gives the comparison results, where we can see that all GCP methods except MPN+E-PN outperform the original network, and all normalization methods improve over the plain COV; among them, our MPN obtains the best result, and outperforms M-Fro, M-$\ell_{2}$ and E-PN by ${\sim}1.3\%$, ${\sim}1.1\%$ and ${\sim}1.3\%$, respectively. Note that GCP with MPN (i.e., MPN-COV) followed by further normalization, either by M-Fro, M-$\ell_{2}$ or E-PN, produces negative gains, so we do not perform any further normalization on our MPN-COV in following experiments.

iSQRT-COV: Number $N$ of Newton-Schulz iterations Here we assess impact of $N$ by employing AlexNet as a backbone model, and using the same hyper-parameters as MPN-COV except the initial learning rate, which is set to $10^{-1.1}$. The top-1 error rate (1-crop prediction) vs. $N$ is illustrated in Fig. 6. With single one iteration, our iSQRT-COV outperforms Plain-COV by $1.3\%$. As iteration number grows, the error rate of iSQRT-COV gradually declines. With 3 iterations, iSQRT-COV is comparable to MPN-COV, having only 0.3% higher error rate, while performing marginally better than MPN-COV between 5 and 7 iterations. The results after $N=7$ show growth of iteration number is not helpful for decreasing error. As larger $N$ incurs higher computational cost, we set $N$ to 5 in the remaining experiments for balancing efficiency and effectiveness. The similar trend was also shown in [24]. An interesting phenomenon is that approximate matrix square root normalization in iSQRT-COV achieves a little gain over the exact one obtained via EIG, which may suggest the fact that matrix square root is not the optimal normalization method for covariance pooling, encouraging development of better normalization methods in the future work.

iSQRT-COV: Significance of post-compensation As discussed in Section 3.4, post-compensation in iSQRT-COV helps to eliminate the side effect resulting from pre-normalization. We verify its significance with ResNet-50 architecture. Table III summarizes impact of different schemes on iSQRT-COV-Net, including simply do nothing (i.e., without post-compensation), Batch Normalization (BN) [60] and our post-compensation scheme. We note that iSQRT-COV-Net fails to converge without post-compensation. Careful observations show that in this case the gradients are very small (on the order of $10^{-5}$), and larger learning rate helps little. Option of BN helps the network converge, but producing about 1% higher top-1 error rate than our post-compensation scheme. The comparison above suggests that our post-compensation is essential for achieving state-of-the-art results with deep CNN architectures, e.g., ResNet-50.

MPN-COV VS. iSQRT-COV in training speed To show acceleration effect of iSQRT-COV, we compare in Fig. 6 training speed between MPN-COV-Net and iSQRT-COV-Net with both 1-GPU and 2-GPU configurations. For 1-GPU configuration, the speed gap vs. batch size between the two methods keeps nearly constant. For 2-GPU configuration, their speed gap becomes more significant when batch size gets larger. As can be seen, the speed of iSQRT-COV-Net continuously grows with increase of batch size while MPN-COV-Net tends to saturate when batch size is larger than 512. Clearly, by avoiding GPU unfriendly EIG or SVD, iSQRT-COV can speed up MPN-COV and make better use of computing power of multiple GPUs.

Compact covariance representations In default setting, our methods output a 32k-dimensional covariance representation. In this paper, the strategies of progressive $1\times 1$ convolutions and group convolution [26] are proposed to further compress covariance representations. We evaluate them using iSQRT-COV with ResNet-50. Table IV summarizes results of our iSQRT-COV and the extra increased parameters (Params.) with respect to ResNet-50 under various compact settings. Specifically, we use a single $1\times 1$ convolution to decrease dimension of convolution features from 2048 to 256, 128 and 64. Then, progressive dimensionality reduction (DR) is performed based on two consecutive $1\times 1$ convolutions, i.e., $2048\rightarrow 512\rightarrow 256$, $2048\rightarrow 512\rightarrow 128$ and $2048\rightarrow 256\rightarrow 64$. We do not employ more $1\times 1$ convolutions as they will bring additional parameters. Meanwhile, size of covariance representation with $2048\rightarrow 512\rightarrow 128$ ($4th$ row) is further reduced to 2K using group convolution of one group, two and four groups, which are indicated by $1G$, $2G$ and $4G$, respectively.

As listed in Table IV, by using a single $1\times 1$ convolution, the recognition error increases about $1.0\%$ and $1.6\%$ when $d$ decreases from 256 to 128 and 64, respectively. The results of progressive DR versions are better than their counterparts based on a single $1\times 1$ convolution, showing the strategy of progressive DR is helpful to preserve classification accuracy for small-size representations. Combining progressive $1\times 1$ convolutions with group convolution, we can obtain $22.4\%$ top-1 error with 2K-dimensional representations, still clearly outperforming the original ResNet-50 with first-order GAP. Above results verify the effectiveness of our compact covariance representations. Our compact strategies can obtain effective small-size representations, very suitable for large scale species classification problem, which will be shown in Section 4.5.

4.2.2 Comparisons with Different Methods

Comparison with other deep global GCP methods We compare different variants of MPN-COV with existing deep global GCP methods using AlexNet architecture, including B-CNN [16], DeepO2P [15] and improved B-CNN [24]. All methods are trained from scratch on ImageNet. We implement B-CNN, DeepO2P and improved B-CNN using public available source code released by the respective authors, and try our best to tune hyper-parameters for them. Note that we use the suggested implementation of improved B-CNN in [24], i.e., FP by SVD and BP by Lyapunov equation. The results of different methods are listed in Table V. We can see that our MPN-COV-Net, G2DeNet and iSQRT-COV-Net obtain similar results, and clearly outperform B-CNN, DeepO2P and improved B-CNN. We owe the gains to matrix power normalization, benefiting robust covariance estimation and proper usage of geometry. The improved B-CNN achieves an unsatisfactory result, which can be regraded as performing element-wise power normalization followed by $\ell_{2}$ normalization behind our MPN-COV. This result suggests that further normalization hurts MPN-COV on large-scale ImageNet classification, which is consistent with the observation in Table II. Additionally, for iSQRT-COV-Net, pre-normalization by trace performs better than Frobenius norm. G2DeNet combining mean vector achieves moderate gains over MPN-COV-Net using larger-size representations, and so we mainly report the results of MPN-COV-Net and iSQRT-COV-Net (trace) in following comparisons.

Why LERM does not work well? DeepO2P exploits matrix logarithm normalization (i.e., LERM [11]) to exploit geometry of covariance, and we claim it is not suitable for deep CNNs. Firstly, the Pow-E metric in MPN-COV improves numerical stability of covariance matrices over LERM. The LERM requires the eigenvalues involved to be strictly positive [11] while the Pow-E metric allows non-negative eigenvalues [18]. For LERM the common method is to add a small positive number $\epsilon$ to eigenvalues for numerical stability. Although $\epsilon$ can be decided by cross-validation, it is difficult to seek the optimal $\epsilon$ suitable for training of deep CNN. For example, [15] suggest $\epsilon=10^{-3}$, which will smooth out eigenvalues less than $10^{-3}$. In contrast, the Pow-E metric do not need such a remedy. Furthermore, from the distribution perspective, the matrix logarithm reverses the order of the significance of eigenvalues, harmful for covariance representations. To make a qualitative analysis, we randomly select 300,000 images from training set of ImageNet, and estimate per-image sample covariances using the outputs of the last convolution layer in AlexNet model. Then, we compute eigenvalues of all covariances using EIG in single-precision floating-point format. The histogram of eigenvalues is shown in Fig. 7 (left), where zero eigenvalues are excluded for better view. Fig. 7 (right) shows the two normalization functions over $[10^{-5},10]$. We can see $\log(\lambda)$ considerably changes the eigenvalue magnitudes, reversing the order of significance of eigenvalues, e.g., a significant eigenvalue $\lambda=50\mapsto\log(\lambda)\approx 3.9$ but an insignificant one $\lambda=10^{-3}\mapsto\log(\lambda)\approx-6.9$. Since significant eigenvalues are generally more important in that they capture the statistics of principal directions along which the feature variances are larger, matrix logarithm will harm covariance pooling.

Comparison with deep local second-order networks In Table VI, we compare our MPN-COV-Net and iSQRT-COV-Net, using ResNet-50 architecture, with deep local second-order networks [46, 47]. Clearly, the two networks integrating local second-order statistics improve over the original one. Our MPN-COV-Net and iSQRT-COV-Net employing matrix square root normalization, are superior to FBN [46] and SORT [47]. These results demonstrate again the effectiveness of our MPN-COV. Similar to the results with AlexNet, iSQRT-COV-Net outperforms MPN-COV-Net by 0.6% in top-1 error.

Comparison with various CNN architectures In the end of this subsection, we evaluate our MPN-COV with different CNN architectures. Since we modify the existing CNN architectures (as described in Section 4.1), we first compare our MPN-COV with the modified CNN architectures. Fig. 8 illustrates the results of different networks with 1-crop prediction, from it we can see that our MPN-COV outperforms consistently the corresponding modified CNN architectures by a clear margin. Then, we combine the proposed MPN-COV with five deep CNN models, including AlexNet [1], VGG-M [64], VGG-VD16 [2], ResNet-50 [4] and ResNet-101 [4]. Here we compare MPN-COV-Net with the original CNN models and the networks with similar architectures. The top-1 error and top-5 error of different networks with 10-crop prediction are given in Fig. 8 and Fig. 8, respectively. The results of compared methods are duplicated from the original papers.

Using AlexNet as backbone model, we compare MPN-COV-Net with the original AlexNet and VGG-F [64]. Our MPN-COV-Net performs much better than both of them. Comparing MPN-COV-Net using VGG-M [64] with the original one, Zeiler & Fergus [65] and OverFeat [66], our MPN-COV-Net shows much better performance than them. Using VGG-VD16 as backbone model, our MPN-COV-Net outperforms the original VGG-VD16 [2] by ${\sim}2.7\%$ in terms of top-1 error, and performs better than GoogleNet [3] and PReLU-net B [67] by ${\sim}1.4\%$ and ${\sim}0.4\%$ in terms of top-5 error, respectively. For ResNet-50 and ResNet-101 architectures, the results show that our MPN-COV-Net performs $1.65\%$ and $2.04\%$ better than the original ones in terms of top-1 error, respectively. Note that our MPN-COV with ResNet-50 and ResNet-101 outperform the original ResNet-101 and ResNet-152 based on first-order GAP, respectively. We also compare iSQRT-COV-Net, in Table VI, with the recently proposed SE-Net [62], CBAM [63] using ResNet-50 and ResNet-101 in 1-crop prediction. Our iSQRT-COV-Net consistently improves SE-Net and CBAM. Compared with DenseNet-201 [5], our iSQRT-COV-Net obtains ${\sim}1.9\%$ gains in top-1 error. Additionally, we make a comprehensive comparison of iSQRT-COV with GAP using ResNets in terms of efficiency and effectiveness, and the corresponding results can be found in Appendix III. These results show our GCP methods can effectively improve CNNs with various architectures with affordable model complexity.

4.3 Large-scale Scene Categorization on Places365

We evaluate our methods on large-scale scene categorization using Places365 dataset [32], which contains about 1.8 million training images and 36,500 validation images collected from 365 scene categories. Following the common settings in [32], we resize all images to $256\times 256$ and randomly crop a $224\times 224$ image patch or its flip for training. The inference is performed with 10-crop prediction, and we report the results on validation set for comparison. Here we compare our iSQRT-COV-Net, using ResNet-50, with five kinds of CNN models, i.e., GoogleNet [3], VGG-VD16 [2], ResNet-50 [4], ResNet-152 [4] and B-CNN with ResNet-50. The results of different methods are given in Table VII, from it we can see that our iSQRT-COV-Net achieves the best results, and outperforms the original ResNet-50 by about $1.2\%$ and $1\%$ in Top-1 and Top-5 errors, respectively. Meanwhile, iSQRT-COV-Net is superior to B-CNN under the same settings. Furthermore, E-PN after our iSQRT-COV block is hurtful in the setting of training from scratch on Places365, which is similar to the results on ImageNet as in Table II. The results on both ImageNet and Places365 verify our proposed methods can significantly improve the representation ability of deep CNNs, achieving much lower classification error.

4.4 Experiments on Small-scale Datasets

In this section, we evaluate our methods on small-scale fine-grained visual recognition and texture classification tasks under the setting of fine-tuning, and assess effect of different fine-tuning strategies.

4.4.1 Fine-grained Visual Recognition

To assess the generalization (transfer) ability of the proposed methods, we pre-train the networks with our GCP on ImageNet and evaluate their performance on fine-grained visual recognition task by fine-tuning. The experiments are conducted on three benchmarks. Among them, Birds-CUB200-2011 [33] is a challenging dataset, including 11,788 images from 200 bird species. FGVC-aircraft [34] is a part of the FGComp 2013 challenge, which consists of 10,000 images across 100 aircraft classes. FGVC-Cars [35] is also presented as a part of the FGComp 2013 challenge, containing 16,185 images from 196 car categories. We employ the fixed training/testing splits provided by the dataset developers, and train or evaluate our networks using neither part annotations nor bounding boxes. For fair comparison, we follow [16] for experimental setting and evaluation protocol. Specifically, we resize the shorter side of input images to 448, and crop center $448\times 448$ patches. We replace 1000-way softmax layer of our pre-trained networks by a k-way softmax layer, where $k$ is number of classes in the corresponding fine-grained dataset, and fine-tune the networks for 50$\sim$100 epochs with a small learning rate $lr$ (e.g., $10^{-2.1}$) for all layers except the last FC layer, which is set to $5\times lr$. The random horizontal flipping is used for data augmentation. After fine-tuning, we perform $\ell_{2}-$normalization on the outputs of our GCP blocks, and feed them to train $k$ one-vs-all linear SVMs with parameter $C=1$. We predict the label of a test image by averaging SVM scores of the image and its horizontal flip.

Table VIII presents classification results of different methods, where our networks significantly improve the original networks under either VGG-VD16 or ResNet-50. When VGG-VD16 is used as backbone model, our proposed methods are superior to both deep CNNs with trainable BoVW (i.e., NetVLAD [49] and NetFV [68]) and deep CNNs with global approximate high-order pooling (i.e., CBP [39], LRBP [40], KP [42], HIHCA [43] and SMSO [45]) by a clear margin. Additionally, our MPN-COV-Net and its variants also outperform other deep GCP networks, i.e., B-CNN [16] and improved B-CNN[24]. With ResNet-50 architecture, iSQRT-COV-Net (8K) respectively outperforms KP by about 2.6%, 3.8% and 0.6% on Birds, Aircrafts and Cars, while iSQRT-COV-Net (32K) further improves accuracy. On all fine-grained datasets, existing methods employing 50-layer ResNet are no better than their counterparts with 16-layer VGG-VD. The reason may be that the last convolution layer of pre-trained ResNet-50 outputs 2048-dimensional features, much higher than 512-dimensional ones of VGG-VD, which are not suitable for existing second- or higher-order pooling methods. Different from all existing methods, our methods perform dimensionality reduction and deem second-order pooling as a component of CNN models, pre-trained on large-scale datasets. Using pre-trained iSQRT-COV-Net with ResNet-101, we establish state-of-the-art results on three fine-grained benchmarks.

4.4.2 Texture Classification

We also transfer our iSQRT-COV-Net to texture classification, where DTD [36] and Indoor67 [37] are employed. Indoor67 has 6,700 images from 67 indoor scene categories, where 80 and 20 images per-category are used for training and test, respectively. DTD consists of 5,640 material images collected from 47 classes, and pre-defined splits in [36] are used for evaluation. We adopt the same experimental settings with [16] for fair comparison. The results of different methods are listed in Table IX, where our iSQRT-COV-Net with VGG-VD16 architecture performs much better than the original model and clearly outperforms other deep second-order pooling networks, i.e., B-CNN [16] and FASON [41]. Meanwhile, our iSQRT-COV-Net is superior to Deep-TEN [69] and MFAFV-Net [52] (i.e., deep BoVW methods) using ResNet-50 and VGG-VD16 as backbone models, respectively. Additionally, iSQRT-COV-Net performs better than SMSO [45] based on both VGG-VD16 and ResNet-50. The results on both FGVC and texture classification demonstrate our GCP methods can significantly improve the generalization (transfer) ability of deep CNNs.

4.4.3 Effect of Different Fine-tuning Strategies.

Additionally, we evaluate how different fine-tuning strategies perform on small-scale datasets using pre-trained models with GAP/GCP on large-scale ImageNet, and whether fine-tuning with E-PN [12, 16] can further improve our proposed methods. Specifically, we employ standard pre-trained networks (i.e., networks with GAP) on ImageNet as backbone models, and replace GAP of the pre-trained models with our iSQRT-COV block while performing fine-tuning on fine-grained and texture benchmarks; for convenience, this method is called iSQRT-COV (PreTr-GAP). As presented in Table X, using both VGG-VD16 and ResNet-50 as backbone models, iSQRT-COV (PreTr-GAP) significantly improves simple fine-tuning of the original GAP-based ones on all benchmarks, indicating that the standard pre-trained networks also can greatly benefit from our GCP (i.e., iSQRT-COV). We also note that iSQRT-COV-Net is superior to iSQRT-COV (PreTr-GAP), which indicates that pre-training the networks with GCP will bring improvement. It also suggests that pre-training the exact models with a large-scale dataset performs better than the pre-training on a surrogate model. Furthermore, we train B-CNN with VGG-VD16 and ResNet-50 as backbone models on ImageNet, and the method is called B-CNN (PreTr-GCP); note that we insert a BN layer after $\ell_{2}$ normalization, otherwise it fails to converge. As compared in Table X, iSQRT-COV clearly outperforms B-CNN under different settings, i.e., iSQRT-COV (PreTr-GAP) and iSQRT-COV-Net are superior to B-CNN and B-CNN (PreTr-GCP), respectively, showing high competitiveness of iSQRT-COV does not solely come from ImageNet pre-training of GCP networks. Finally, we pre-train the networks with iSQRT-COV and perform fine-tuning with E-PN, and this method is indicated by iSQRT-COV + E-PN. The results of iSQRT-COV-Net vs. iSQRT-COV + E-PN show that E-PN has little effect on our iSQRT-COV under the setting of fine-tuning.

4.5 iNaturalist Challenge 2018

As part of the FGVC5 workshop at CVPR 2018, iNaturalist Challenge 2018 is a large-scale species classification competition, which contains over 8,000 species with 437,513 training, 24,426 validation and 149,394 test images. This dataset suffers from many visually similar species and high class imbalance, which are extremely difficult for accurate classification without expert knowledge. We use ResNet-152 pre-trained on ImageNet-11K as backbone model, and perform fine-tuning on training and validation sets. The top-3 error is reported on test server for comparison. To efficiently classify a large number of species, we use fast MPN-COV (i.e., iSQRT-COV with 3 iterations), and use a $1\times 1$ convolution with $160$ channels before GCP, leading to $\sim$12K dimensional representations, connected to 8142-way softmax classifier. Table XI gives the results of compared methods, where iSQRT-COV-Net distinctly outperforms the original ResNet-152 in configurations of both different input image sizes and ensemble models. By fusing iSQRT-COV-Net and ResNet-152 models, we positioned the first place, and outperforms the runner-up by ${\sim}1.2\%$. Note that only ensemble of iSQRT-COV-Net models is still better than the runner-up by ${\sim}0.8\%$, showing effectiveness of our method.

5 Conclusion

This paper proposed a global Matrix Power Normalized COVariance (MPN-COV) Pooling methodology for improving the representation and generalization abilities of deep CNNs. Our matrix power normalization can not only estimate robustly covariance matrices but also can properly exploit Riemannian geometry of covariances. To further improve MPN-COV, we proposed a global Gaussian embedding method to integrate additional first-order statistical information, developed an iterative matrix square root normalization method for speeding up training of our covariance pooling networks, and studied a compact covariance representation strategy to reduce model complexity. As a result, our matrix power normalization methodology can improve significantly existing global covariance pooling methods, while enjoying fast training speed and affordable model complexity; meanwhile, it has solid statistical and geometrical foundations. The comprehensive experiments on both large-scale and small-scale benchmarks verified our MPN-COV and its improved solutions achieve a clear improvement with diversified deep CNN architectures, while obtaining state-of-the-art classification performance. In view of the effectiveness and architecture-independence of the proposed methods, they have the potential to be a commodity component in the existing deep CNN models, and they can motivate interests in future exploration of higher-order information in design of novel CNN architectures. In future, we will apply the proposed methods to other vision tasks, such as object detection [70] and semantic segmentation [71].

Acknowledgments

The work was supported by the National Natural Science Foundation of China (Grant No. 61971086, 61806140, U19A2073, 61471082, 61671182 and 61732011). Q. Wang was supported by National Postdoctoral Program for Innovative Talents. The work was done while Q. Wang was a PhD student at the School of Information and Communication Engineering, Dalian University of Technology, China.

Appendix I: Derivations for Backpropagation of iSQRT-COV Block

Though autograd toolkits provided by some deep learning frameworks can accomplish backpropagation of iSQRT-COV block automatically, their involved BP is still in a black box. Meanwhile, autograd toolkits sometimes bring uncertainty, e.g., autograd toolkit of PyTorch 0.3.0 or below cannot compute gradients of iSQRT-COV correctly. To make iSQRT-COV be self-contained and enable its implementation to be accessible when autograd toolkit is unavailable (e.g., the well-known Caffe and early versions of MatConvNet), we derive the gradients associated with the structured layers based on matrix backpropagation methodology [55]. Below we take pre-normalization by trace as an example, deriving the corresponding gradients.

BP of Post-compensation Given $\frac{\partial l}{\partial\mathbf{Z}}$ where $l$ is the loss function, the chain rule is of the form \mathrm{tr}\big{(}\big{(}\frac{\partial l}{\partial\mathbf{Z}}\big{)}^{T}\mathrm{d}\mathbf{Z}\big{)}=\mathrm{tr}\big{(}\big{(}\frac{\partial l}{\partial\mathbf{Y}_{N}}\big{)}^{T}\mathrm{d}\mathbf{Y}_{N}+\big{(}\frac{\partial l}{\partial\boldsymbol{\Sigma}}\big{)}^{T}\mathrm{d}\boldsymbol{\Sigma}\big{)}. After some manipulations, we have

[TABLE]

BP of Newton-Schulz Iteration Then we compute the partial derivatives with respect to $\frac{\partial l}{\partial\mathbf{Y}_{k}}$ and $\frac{\partial l}{\partial\mathbf{P}_{k}}$, $k=N-1,\ldots,1$, given $\frac{\partial l}{\partial\mathbf{Y}_{N}}$ obtained by Eq. (Appendix I: Derivations for Backpropagation of iSQRT-COV Block) and $\frac{\partial l}{\partial\mathbf{P}_{N}}=0$. As the matrix $\boldsymbol{\Sigma}$ is symmetric, it is easy to see from equation of Newton-Schulz iteration (i.e., Eq. (3.4)) that $\mathbf{Y}_{k}$ and $\mathbf{P}_{k}$ are both symmetric. According to the chain rules (omitted hereafter for simplicity) of matrix backpropagation and after some manipulations, for $k=N,\ldots,2$, we can derive

[TABLE]

The final step of this layer is concerned with the partial derivative with respect to $\frac{\partial l}{\partial\mathbf{A}}$, which is given by

[TABLE]

BP of Pre-normalization Note that here we need to combine the gradient of the loss function $l$ with respect to $\boldsymbol{\Sigma}$, backpropagated from the post-compensation layer. As such, by referring to equation of pre-normalization layer, we make similar derivations as before and obtain

[TABLE]

If we adopt pre-normalization by Frobenius norm, the gradients associated with post-compensation become

[TABLE]

and that with respect to pre-normalization is

[TABLE]

while the backward gradients of Newton-Schulz iteration (Appendix I: Derivations for Backpropagation of iSQRT-COV Block) keep unchanged. Based on above derivations, we can achieve the backpropagation of iSQRT-COV block without autograd toolkit.

Appendix II: Comparison of Different Gaussian Embedding Methods for G2DeNet

To insert a global Gaussian distribution into deep CNNs, we introduce a Gaussian embedding method studied in [23] to identify a Gaussian $\mathcal{N}(\boldsymbol{\mu},\boldsymbol{\Sigma})$ as a square root SPD matrix (i.e., $\left[\begin{smallmatrix}\boldsymbol{\Sigma}+\boldsymbol{\mu}\boldsymbol{\mu}^{T}&\boldsymbol{\mu}\\ \boldsymbol{\mu}^{T}&1\end{smallmatrix}\right]^{\frac{1}{2}}$). Although many Gaussian embedding methods have been studied, there exist clear differences between ours and them. Specifically, Nakayama et al. [59] embed Gaussians in a flat manifold by taking an affine coordinate system. In [58], Gaussian is mapped to a unique positive definite lower triangular affine transform (PDLTAT) matrix (i.e., $\left[\begin{matrix}\mathbf{L}&\boldsymbol{\mu}\\ \mathbf{0}^{T}&1\end{matrix}\right]$ and $\mathbf{L}\mathbf{L}^{T}$ is the Cholesky decomposition of $\boldsymbol{\Sigma}$), whose space forms an affine group. Such Gaussian embedding is not suitable for backpropagation due to involvement of Cholesky decomposition of covariance matrix $\boldsymbol{\Sigma}$. The methods in Calvo et al. [56] and Lovri’c et al. [57] respectively embed the space of Gaussian in the Siegel group and the Riemannian symmetric space, identifying a Gaussian as a unique SPD matrix, i.e., $\begin{bmatrix}\boldsymbol{\Sigma}+\boldsymbol{\mu}\boldsymbol{\mu}^{T}&\boldsymbol{\mu}\\ \boldsymbol{\mu}^{T}&1\end{bmatrix}$ or $C\ast\begin{bmatrix}\boldsymbol{\Sigma}+\boldsymbol{\mu}\boldsymbol{\mu}^{T}&\boldsymbol{\mu}\\ \boldsymbol{\mu}^{T}&1\end{bmatrix}$ where $C$ is a scalar. After that, [72, 73] employ Log-Euclidean Riemannian metric (LERM) [11] on the resulting embedding matrix in [56, 57], which can be presented as \log\bigg{(}\begin{bmatrix}\boldsymbol{\Sigma}+\boldsymbol{\mu}\boldsymbol{\mu}^{T}&\boldsymbol{\mu}\\ \boldsymbol{\mu}^{T}&1\end{bmatrix}\bigg{)}. Different from these methods, our introduced embedding method considers both geometric and algebraic (i.e., Lie group) structures of Gaussian, while achieving better performances. In particular, we compare our G2DeNet with various Gaussian embedding methods on Birds-CUB200-2011 dataset using VGG-VG16 as backbone model. The results are given in Table XII, from it we can see that our introduced embedding method achieves the best performance, outperforming other competing methods by $3\%\sim 3.6\%$.

Appendix III: Comparison of iSQRT-COV and GAP in terms of efficiency and effectiveness using ResNets

In this section, we compare global average pooling (GAP) and our proposed global covariance pooling (GCP) on large-scale ImageNet using ResNet-18, ResNet-34, ResNet-50 and ResNet-101 as backbone models. The evaluation metrics include network parameters, floating point operations per second (FLOPs), inference (FP computation) time per image, memory cost for testing one image and Top-1/Top-5 classification errors, given input images of size 224$\times$224. For our GCP, we use the iSQRT-COV-Net by setting covariance representations to 32k and 8k, respectively. All models are trained with the same experimental settings and run on a workstation equipped with four Titan Xp GPUs, two Intel(R) Xeon Silver 4112 CPUs @ 2.60GHz, 64G RAM and 480 GB INTEL SSD. From the comparison results in Table XIII, we have three observations as follows. (1) Both our iSQRT-COV (32k) and iSQRT-COV (8k) significantly outperform the original first-order GAP under the same CNN architectures. In particular, iSQRT-COV (8k) achieves clear gains with affordable computational cost. Specifically, it introduces extra $\sim$7M parameters, $\sim$0.2ms inference time and $\sim$100M memory cost, but decreases about 5.3%, 3.5%, 2.4% and 2.3% Top-1 errors over the original GAP-based ResNet-18, ResNet-34, ResNet-50 and ResNet-101, respectively. (2) When we compare different models that have matching performance, we can see that iSQRT-COV (8k) models show much lower complexity than the corresponding GAP models. For example, ResNet34 + iSQRT-COV (8k) and ResNet50 + iSQRT-COV (8k) are comparable or slightly superior to ResNet101 + GAP and ResNet152 + GAP in performance, respectively, whereas the former ones have much lower complexity than the latter ones. (3) When different models that have similar space and time complexity are compared, it can be seen that iSQRT-COV (8k) models clearly outperform the corresponding GAP models. For example, ResNet18 + iSQRT-COV (8k), ResNet50 + iSQRT-COV (8k) and ResNet101 + iSQRT-COV (8k) have similar (or lower) complexity with ResNet34 + GAP, ResNet101 + GAP and ResNet152 + GAP, but drop 1.8%, 1.3% and 1.7% in terms of Top-1 error, respectively. Above results demonstrate that our iSQRT-COV has better ability to balance efficiency and effectiveness.

Appendix IV: Comparison between SVM and Softmax classifiers on Fine-grained Visual Classification Task

For the purpose of fair comparison with existing works, we in fact adopt exactly the same experimental settings with [16, 39, 42, 43], where an SVM is trained for classification after fine-tuning models. To assess the effect of classifiers on our method, we provide the results with softmax classifier used at training time, as presented in Table XIV. We can see that SVM and softmax classifiers achieve very similar results for our iSQRT-COV-Net, except iSQRT-COV-Net with ResNet-50 on Aircrafts, where softmax classifier is a little better than SVM classifier. The experiments suggest that, for our iSQRT-COV-Net, overall training SVM classifier after fine-tuning with softmax classifier brings negligible benefits.

Appendix V: Comparison of iSQRT-COV and Improved B-CNN in terms of Memory, Speed and Accuracy

Here, we compare iSQRT-COV with Improved B-CNN [24] in term of memory cost, training speed and classification accuracy. Our induces memory overhead as the intermediate tensors in the Newton-Schulz iteration have to be cached for usage in the backpropagation; by contrast, the methods proposed in improved B-CNN [24] do not involve this memory overhead. We conduct additional experiments to evaluate speed and memory cost of different methods. The improved B-CNN proposed several different fine-tuning methods, including FP and BP both dependent upon SVD, FP by Denman-Beavers iteration or Newton-Schulz (NS) iteration while BP by SVD$+$Lyapunov (Lyap) equation. We compare with the method that uses NS iteration for FP and SVD$+$Lyap for BP, which is more efficient among all the proposed methods and allow accurate gradients computation. The authors implement, in their released code, a new method that was not described in their paper, in which gradients can be computed by iterative Lyap solver, avoiding time consuming SVD; we also compare with this method. We make experiments on Birds-CUB200-2011 dataset following the experimental settings of improved B-CNN [24]. We employ VGG-VD16 as backbone model and set batch-size to 24. To implement improved B-CNN, we use the source code released by the authors. We employ the function ’sqrt_newton_schulz’ for FP by NS iteration, the function ’sqrt_svd_lyap’ for BP by SVD$+$Lyap, and the function ’lyap_newton_schulz’ for BP by iterative Lyap solver. The experiments run on a PC equipped with a 4-core Intel i7-4790k@4.0GHz CPU, 32G RAM, 512GB Samsung PRO SSD and a single GTX 1080Ti. The results are presented in Table XV, where iSQRT-COV brings extra memory (i.e., $+\Delta$) as it needs to cache intermediate tensors. Nevertheless, the percentage of memory overhead with respect to that of the whole networks (ALL) is negligible ($<$1%). The training speed of iSQRT-COV is 1.8x faster than improved B-CNN based on NS iteration and SVD+lyap. Note that, as shown in project page of the authors, iterative Lyapunov solver is an approximation algorithm, leading to larger error. Therefore, it requires more iterations for achieving satisfying results. When number of iterations is 5, improved B-CNN based on NS iteration and Lyap iteration has comparable training speed with iSQRT-COV, but iSQRT-COV has higher classification accuracy. More iterations will bring some improvement for improved B-CNN based on NS iteration and Lyap iteration, but it slows down training speed. Overall, iSQRT-COV can make a good balance for memory cost, training speed and classification accuracy.