Low-Rank plus Sparse Decomposition of Covariance Matrices using Neural Network Parametrization

TL;DR

This paper proposes a neural network-based method for decomposing positive semidefinite matrices into low-rank and sparse components, with theoretical convergence guarantees and applications in portfolio optimization.

Contribution

It introduces a neural network parametrization for the low-rank component and analyzes its convergence rate, advancing matrix decomposition techniques.

Findings

Convergence rate grows polynomially with input and network dimensions.

Method effectively decomposes covariance matrices in portfolio optimization.

Provides theoretical analysis of gradient descent convergence.

Abstract

This paper revisits the problem of decomposing a positive semidefinite matrix as a sum of a matrix with a given rank plus a sparse matrix. An immediate application can be found in portfolio optimization, when the matrix to be decomposed is the covariance between the different assets in the portfolio. Our approach consists in representing the low-rank part of the solution as the product $MM^{T}$, where $M$ is a rectangular matrix of appropriate size, parametrized by the coefficients of a deep neural network. We then use a gradient descent algorithm to minimize an appropriate loss function over the parameters of the network. We deduce its convergence rate to a local optimum from the Lipschitz smoothness of our loss function. We show that the rate of convergence grows polynomially in the dimensions of the input, output, and the size of each of the hidden layers.

Tables

Table 1. TABLE I : Choice of activation function σ 𝜎 \sigma .

Name	Definition	${𝝈}_{\mathrm{𝐦𝐚𝐱}}^{\prime }$	${𝝈}_{\mathrm{𝐦𝐚𝐱}}^{\prime \prime }$
Logistic	$\frac{1}{1+e^{-x}}$	0.25	$\frac{1}{6\sqrt{3}}$
Hyperbolic tangent	$\frac{e^{2x}-1}{e^{2x}+1}$	1	$\frac{4}{3\sqrt{3}}$
(Scaled) $\arctan$	$\frac{2}{\pi }{\tan }^{-1}\left(x\right)$	$\frac{2}{\pi }$	$\frac{3\sqrt{3}}{4\pi }$

Table 2. TABLE II: Comparison between our method (DNN) and several other algorithms: IALM [ 13 ] , Non-convex RCPA (NC-RCPA) [ 12 ] , RPCA-GD [ 27 ] , and FPCP [ 21 ] . The input parameters are the dimension n 𝑛 n of the matrix Σ Σ \Sigma the given rank k 0 subscript 𝑘 0 k_{0} of L 0 subscript 𝐿 0 L_{0} , the given sparsity s 0 subscript 𝑠 0 s_{0} of S 0 subscript 𝑆 0 S_{0} , and the chosen forced rank k 𝑘 k of L 𝐿 L . We report the estimated rank r ​ ( L ) 𝑟 𝐿 r(L) of the output matrix L 𝐿 L , the estimated sparsity s ​ ( S ) 𝑠 𝑆 s(S) of the output S 𝑆 S , and their respective relative errors.

$n$	$k_0$	$k$	$s_0$	Algorithm	$r\left(L\right)$	$s\left(S\right)$	$\text{rel.error}\left(S\right)$	$\text{rel.error}\left(L\right)$
100	10	5	0.60	IALM	11.00 (0.45)	0.17 (0.02)	0.98 (0.00)	0.93 (0.01)
				NC-RCPA	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.95 (0.01)
				RPCA-GD	5.00 (0.00)	0.01 (0.00)	0.99 (0.00)	0.94 (0.01)
				FPCP	5.00 (0.00)	0.04 (0.00)	1.01 (0.01)	0.97 (0.01)
				DNN	5.00 (0.00)	0.01 (0.00)	0.36 (0.02)	0.53 (0.02)
			0.95	IALM	11.30 (0.64)	0.16 (0.01)	2.21 (0.16)	0.28 (0.02)
				NC-RCPA	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.13 (0.00)
				RPCA-GD	5.00 (0.00)	0.01 (0.00)	4.18 (0.15)	0.53 (0.02)
				FPCP	5.00 (0.00)	0.04 (0.00)	6.41 (0.21)	0.84 (0.01)
				DNN	5.00 (0.00)	0.01 (0.00)	4.35 (0.18)	0.54 (0.02)
		10	0.60	IALM	11.00 (0.45)	0.17 (0.02)	0.98 (0.00)	0.93 (0.01)
				NC-RCPA	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.95 (0.01)
				RPCA-GD	10.00 (0.00)	0.03 (0.01)	0.99 (0.00)	0.94 (0.01)
				FPCP	9.00 (0.00)	0.07 (0.01)	0.99 (0.01)	0.95 (0.01)
				DNN	10.00 (0.00)	0.06 (0.02)	0.03 (0.01)	0.04 (0.01)
			0.95	IALM	11.30 (0.64)	0.16 (0.01)	2.21 (0.16)	0.28 (0.02)
				NC-RCPA	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.13 (0.00)
				RPCA-GD	10.00 (0.00)	0.09 (0.18)	1.87 (0.38)	0.24 (0.05)
				FPCP	9.00 (0.00)	0.07 (0.00)	5.64 (0.20)	0.74 (0.01)
				DNN	10.00 (0.00)	0.16 (0.05)	0.14 (0.05)	0.02 (0.01)
		20	0.60	IALM	11.00 (0.45)	0.17 (0.02)	0.98 (0.00)	0.93 (0.01)
				NC-RCPA	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.95 (0.01)
				RPCA-GD	20.00 (0.00)	0.06 (0.01)	0.97 (0.00)	0.92 (0.01)
				FPCP	12.00 (0.00)	0.19 (0.02)	0.97 (0.00)	0.93 (0.01)
				DNN	20.00 (0.00)	0.06 (0.01)	0.05 (0.01)	0.08 (0.01)
			0.95	IALM	11.30 (0.64)	0.16 (0.01)	2.21 (0.16)	0.28 (0.02)
				NC-RCPA	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.13 (0.00)
				RPCA-GD	20.00 (0.00)	0.06 (0.01)	2.03 (0.18)	0.26 (0.02)
				FPCP	12.00 (0.00)	0.19 (0.01)	3.99 (0.19)	0.53 (0.01)
				DNN	20.00 (0.00)	0.17 (0.04)	0.32 (0.08)	0.04 (0.01)
200	5	5	0.95	IALM	5.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.27 (0.01)
				NC-RCPA	5.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.27 (0.01)
				RPCA-GD	5.00 (0.00)	0.99 (0.00)	1.00 (0.00)	0.27 (0.01)
				FPCP	5.00 (0.00)	0.06 (0.00)	2.90 (0.14)	0.79 (0.01)
				DNN	5.00 (0.00)	0.18 (0.02)	0.08 (0.01)	0.02 (0.00)
400	10	10	0.95	IALM	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.27 (0.01)
				NC-RCPA	10.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.27 (0.01)
				RPCA-GD	10.00 (0.00)	0.99 (0.00)	1.00 (0.00)	0.27 (0.01)
				FPCP	7.50 (0.50)	0.02 (0.00)	3.13 (0.07)	0.87 (0.00)
				DNN	10.00 (0.00)	0.12 (0.05)	0.08 (0.03)	0.02 (0.01)
800	20	20	0.95	IALM	20.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.28 (0.00)
				NC-RCPA	20.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.28 (0.00)
				RPCA-GD	20.00 (0.00)	1.00 (0.00)	1.00 (0.00)	0.28 (0.00)
				FPCP	9.00 (0.00)	0.01 (0.00)	3.30 (0.03)	0.92 (0.00)
				DNN	18.10 (5.70)	0.04 (0.02)	1.16 (1.59)	0.32 (0.44)

Equations

$$\underset{L,S}{\min }||L|{|}_{\ast }+\delta ||S|{|}_{1}s.t.\Sigma =L+S.$$

$$\underset{L,S,Y}{\min }||L|{|}_{\ast }+\delta ||S|{|}_1+⟨Y,\Sigma -L-S⟩+\frac{\mu }{2}||\Sigma -L-S|{|}_F^2,$$

$$L_{t+1}$$

$$S_{t+1}$$

$$Y_{t+1}$$

$${\mu }_{t+1}$$

$$\underset{L,S}{\min }\delta ||S|{|}_1+\frac{1}{2}||\Sigma -L-S|{|}_F^{2}s.t.\mathrm{rank}(L)=k.$$

$$L_{t+1}$$

$$S_{t+1}$$

$$h:{\mathbb{S}}^n$$

$$\Sigma$$

$$f_u^{A^{(u)},b^{(u)}}(v):=A^{(u)}v+b^{(u)}.$$

$$\begin{array}{r}{f}_m^{A^{(m)},b^{(m)}}\circ {\sigma }^{(m)}\circ f_{m-1}^{A^{(m-1)},b^{(m-1)}}\circ \dots \circ {\sigma }^{(1)}\circ f_0^{A^{(0)},b^{(0)}}(x).\end{array}$$

$$\varphi_{\mathrm{obj}}(\Theta):=\left\|g\big{(}\mathcal{N}^{\Theta}_{m}\left(h(\Sigma)\right)\!\big{)}\,g\big{(}\mathcal{N}^{\Theta}_{m}\left(h(\Sigma)\right)\!\big{)}^{T}-\Sigma\right\|_{1}.$$

$${\varphi}(\Theta):=\sum_{i,j=1}^{n}\mu\left(\left[g\big{(}\mathcal{N}^{\Theta}_{m}\left(h(\Sigma)\right)\!\big{)}\,g\big{(}\mathcal{N}^{\Theta}_{m}\left(h(\Sigma)\right)\!\big{)}^{T}-\Sigma\right]_{i,j}\right),$$

$$\mu (t):=\{\begin{array}{ll}\frac{t^2}{2\epsilon }+\frac{\epsilon }{2}& \text{if }|t|\le \epsilon \\ |t|& \text{if }|t|>\epsilon ,\end{array}$$

$$\begin{split}\begin{array}[]{l}\textrm{Fix}\;\Theta_{0}\\ \textrm{For}\;j\geq 0:\\ \quad\quad\textrm{Compute }\nabla{\varphi}(\Theta_{j})\\ \quad\quad\textrm{Determine a step-size }h_{j}>0\\ \quad\quad\text{Set}\;\Theta_{j+1}=\Theta_{j}-h_{j}\nabla{\varphi}(\Theta_{j}).\end{array}\end{split}$$

$$\Vert \left(X^1,\dots ,X^{\gamma }\right)\Vert :={\left(||X^1|{|}_F^2+\cdots +||X^{\gamma }|{|}_F^2\right)}^{1/2}.$$

$$\underset{j\in {\mathbb{N}}_0}{\sup }\Vert {\Theta }_j\Vert \le D.$$

$$\displaystyle L^{2}_{m}=\mathcal{C}_{m}\max\Big{\{}kn^{3}D^{4}\ell_{\max},$$

$$\displaystyle\max\{nD^{2}||h(\Sigma)||^{2},nD^{2}\ell_{\max}\Big{\}},$$

$$\begin{split}\hat{L}^{2}_{1}=\mathcal{C}_{1}n^{2}\ell_{1}\max\!\bigg{\{}&D^{4}||h(\Sigma)||^{4},\ell_{1}D^{8}||h(\Sigma)||^{4},\\ &\ell_{1}^{2}D^{6}||h(\Sigma)||^{2},nD^{2}||h(\Sigma)||^{2},nk\ell_{1},\ell_{1}^{3}D^{4}\bigg{\}}.\end{split}$$

$$\phi ({\Theta }_j)-\phi ({\Theta }_{j+1})\ge \frac{K}{L_m}\Vert \nabla \phi ({\Theta }_j){\Vert }^2,$$

$$\min_{0\leq j\leq N}\|\nabla\varphi(\Theta_{j})\|\leq\tfrac{1}{\sqrt{N+1}}\Big{[}\tfrac{L_{m}}{K}\big{(}\varphi(\Theta_{0})-\varphi^{*}\big{)}\Big{]}^{1/2},$$

$$N+1\geq\tfrac{L_{m}}{K\varepsilon^{2}}\big{(}\varphi(\Theta_{0})-\varphi^{*}\big{)}\ \ \Longrightarrow\ \ \min_{0\leq j\leq N}\|\nabla\varphi(\Theta_{j})\|\leq\varepsilon.$$

$$\begin{split}\alpha\big{\langle}\nabla\varphi(\Theta_{j}),\Theta_{j}-\Theta_{j+1}\big{\rangle}&\leq\varphi(\Theta_{j})-\varphi(\Theta_{j+1})\\ \beta\big{\langle}\nabla\varphi(\Theta_{j}),\Theta_{j}-\Theta_{j+1}\big{\rangle}&\geq\varphi(\Theta_{j})-\varphi(\Theta_{j+1}).\end{split}$$

$$L_m\le O\left(k^{1/2}{n}^{3/2}{D}^{2m+2}{\ell }_{\max }||h(\Sigma )|{|}^2\right).$$

$$\underset{0\le j\le N}{\min }\Vert \nabla \phi ({\Theta }_j)\Vert \le \frac{1}{\sqrt{N+1}}{\left(\frac{L_m}{K}\phi ({\Theta }_0)\right)}^{1/2}.$$

$$f_1^{A,b}:{\mathbb{R}}^r\to {\mathbb{R}}^{\ell },w\mapsto f_1^{A,b}(w)=Aw+b.$$

$$f_2^{C,d}:{\mathbb{R}}^{\ell }\to {\mathbb{R}}^{nk}v\mapsto f_2^{C,d}(v)=Cv+d.$$

$$\tilde{\phi }(\Theta )=\sum _{i,j=1}^n\mu \left({\left[g\left({\mathcal{N}}^{\Theta }\left(h(\Sigma )\right)\right)g{\left({\mathcal{N}}^{\Theta }\left(h(\Sigma )\right)\right)}^T-\Sigma \right]}_{i,j}\right),$$

$$\begin{array}{rl}{\omega }_{i,j}& :={\left[g\left({\mathcal{N}}^{\Theta }\left(h(\Sigma )\right)\right)g{\left({\mathcal{N}}^{\Theta }\left(h(\Sigma )\right)\right)}^T-\Sigma \right]}_{i,j},\\ X& :={\mathcal{N}}^{\Theta }\left(h(\Sigma )\right)\in {\mathbb{R}}^{nk},\\ Y& :=\tilde{\sigma }\circ f_1^{A,b}(h(\Sigma ))\in {\mathbb{R}}^{\ell },\\ Z& :=f_1^{A,b}(h(\Sigma ))\in {\mathbb{R}}^{\ell }.\end{array}$$

$$\begin{array}{rl}\frac{\partial \tilde{\phi }(\Theta )}{\partial b_{\iota }}=\sum _{i,j=1}^n& {\mu }^{\prime }({\omega }_{i,j}){\sigma }^{\prime }(Z_{\iota })\times \\ & \left(\sum _{s=1}^k\left(C_{(i-1)k+s,\iota }{X}_{(j-1)k+s}+X_{(i-1)k+s}{C}_{(j-1)k+s,\iota }\right)\right).\end{array}$$

Full text

I Introduction

We present a new approach to decompose a possibly large covariance matrix $\Sigma$ into the sum of a positive semidefinite low-rank matrix $L$ plus a sparse matrix $S$. Our approach consists in fixing an (upper bound for the) rank $k$ of $L$ by defining $L:=MM^{T}$ for a suitable $M\in\mathbb{R}^{n\times k}$, where one parametrizes $M$ using a deep neural network, whose coefficients are minimized using a gradient descent method.

Albeit our method can be used in any context where such a problem occurs, our primary application of interest is rooted in Finance. When studying the correlation matrix between the returns of financial assets, it is important for the design of a well-diversified portfolio to identify groups of heavily correlated assets, or more generally, to identify a few ad-hoc features, or economic factors, that describe some dependencies between these assets. To this effect, the most natural tool is to determine the few first dominant eigenspaces of the correlation matrix and to interpret them as the main features driving the behavior of the portfolio. This procedure, generally termed Principal Component Analysis (PCA), is widely used. However, PCA does not ensure any sparsity between the original matrix $\Sigma$ and its approximation $A$. As it turns out, many coefficients of $\Sigma-A$ can be relatively large with respect to the others; these indicate pairs of assets that present an ignored large correlation between themselves, beyond the dominant features revealed by PCA. Following [24], to reveal this extra structure present in $\Sigma$, we decompose it into the sum of a low-rank matrix $L$, which describes the dominant economic factors of a portfolio, plus a sparse matrix $S$, to identify hidden large correlations between assets pairs. The number of those economic factors is set according to the investor’s views on the market, and coincides with the rank of $L$. The sparse part $S$ can be seen as a list of economic abnormalities, which can be exploited by the investor.

Beyond covariance matrices, this decomposition is a procedure abundantly used in image and video processing for compression and interpretative purposes [3], but also in latent variable model selection [6], in latent semantic indexing [16], in graphical model selection [2], in graphs [23], and in gene expression [14], among others. A rich collection of algorithms exist to compute such decomposition, see [7, 9, 13, 25, 26] to cite but a few, most of which are reviewed in [4] and implemented in the Matlab LRS library [5]. However, our method can only address the decomposition of a symmetric positive semidefinite matrix, as it uses explicitly and takes full advantage of this very particular structure.

The Principal Component Pursuit (PCP) reformulation of this problem has been proposed and studied by [6, 7, 13] as a robust alternative to PCA, and generated a number of efficient algorithms. For a given $\delta>0$, the PCP problem is formulated as

[TABLE]

where $\Sigma\in\mathbb{R}^{n\times n}$ is the observed matrix, $||L||_{*}$ is the nuclear norm of matrix $L$ (i.e. the sum of all singular values of $L$) and $||S||_{1}$ is the $l^{1}$-norm of matrix $S$. To solve (1), an approach consists in constructing its Augmented Lagrange Multiplier (ALM) [10]. By incorporating the constraints of (1) into the objective multiplied by their Lagrange multiplier $Y\in\mathbb{R}^{n\times n}$, the ALM is

[TABLE]

which coincides with the original problem when $\mu\to\infty$. We denote by $\mathcal{L}_{\mu}(L,S,Y)$ the above objective function. In [13], it is solved with an alternating direction method:

[TABLE]

for some $\rho>1$. The resulting method is called Inexact ALM (IALM); in the Exact ALM, $\mathcal{L}_{\mu}(L,S,Y_{t})$ is minimized on $L$ and $S$ simultaneously, a considerably more difficult task. The problem (4) can be solved explicitly at modest cost. In contrast, (3) requires an expensive Singular Value Decomposition (SVD). In [12], the authors replace the nuclear norm in (2) by the non-convex function $||L||_{\gamma}:=\sum_{i}(1+\gamma)\sigma_{i}(L)/(\gamma+\sigma_{i}(L))$ that interpolates between $\operatorname{rank}(L)$ and $||X||_{*}$ as $\gamma$ goes from [math] to $1$, in order to depart from the convex PCP approximation of the original problem. Then they apply the same alternating direction method to the resulting function. This method is referred to as Non-Convex RCPA or NC-RCPA. In [21], the authors rather solve a variant of (1) by incorporating the constraint into the objective, removing the costly nuclear norm term, and imposing a rank constraint on $L$:

[TABLE]

We denote by $\widehat{\mathcal{L}}(L,S)$ the above objective function. Using also an alternating direction strategy, the authors have devised the Fast PCP (FPCP) method as

[TABLE]

The problem (7) is easy to solve as for (4). In contrast to (3), the sub-problem (6) can be solved by applying a faster partial SVD to $\Sigma-S_{t}$, with the only necessity of computing the $k$ first singular values and their associated eigenvectors. These authors have further improved their algorithm in [22] with an accelerated variant of the partial SVD algorithm. Their methods are considered as state-of-the-art in the field [4] and their solution is of comparable quality to the one of (1). An alternative approach, designated here as (RPCA-GD) to solve (5) was proposed in [27], where, as in our setting, the rank constraint is enforced by setting $L:=MM^{T}$ for $M\in\mathbb{R}^{n\times k}$. Then a projected gradient method is used to solve (6) in $M$. In order to guarantee that $S$ has a prescribed sparsity, they use an ad-hoc projector on an appropriate space of sparse matrices.

The solution to the PCP problem (1) depends on the hyperparameter $\delta$, from which we cannot infer the value of the rank $k$ of the resulting optimal $L$ with more accuracy than in the classical results of Candès et al. in Theorem 1.1 of [6]. In view of the particular financial application we have in mind, we would prefer a method for which we can chose the desired $k$ in advance. In both the IALM and the Non-Convex RPCA methods, neither the rank of $L$ nor its expressivity — that is, the portion of the spectrum of $\Sigma$ covered by the low-rank matrix — can be chosen in advance. In RPCA-GD, the rank is chosen in advance, but the sparsity of $S$ is set in advance by inspection of the given matrix $\Sigma$, a limitation that we would particularly like to avoid in our application. From this point of view, FPCP seems the most appropriate method. In our approach, we must first select a rank for $L$, based e.g. on a prior spectral decomposition of $\Sigma$ or based on exogenous considerations. We then apply a gradient descent method with a well-chosen loss function, using Tensorflow [1] or Pytorch [20].

In Section II, we introduce the construction of our low-rank matrix $L=MM^{T}$, where $M$, in contrast with the RPCA-GD method, is parametrized by the coefficients of a multi-layered neural network. A potential advantage of this parametrization has been pointed in [15], albeit in a different context than ours. We also describe the corresponding loss function that we seek to minimize. Moreover, we analyze the regularity properties of the objective function leading to an estimate of the convergence rate of a standard gradient descent method to a stationary point of the method; see Theorem II.1. We show that the convergence rate of our algorithm grows polynomially in the dimension of each layer. In Section III, we conduct a series of experiments first on artificially generated data, that is matrices $\Sigma$ with a given decomposition $L+S$, to assess the accuracy of our method and to compare it with the algorithms presented in this section. We apply also our algorithm to real data sets: the correlation matrix of the stocks from the S&P500 index and an estimate of the correlation matrix between real estate market indices of 44 countries. We show that our method achieves a higher accuracy than the other algorithms. Moreover, by its construction of $L:=MM^{T}$, we can guarantee the positive semidefiniteness of $L$, although empirical covariance matrices tend to not satisfy this property. We refer to [11] for a detailed discussion of this issue. By contrast, most algorithms do not ensure that $L$ is kept positive semidefinite, which forces them to correct their output at the expense of their accuracy. The proof of Theorem II.1 is provided in Section IV.

II Neural network parametrized optimization and its convergence rate

Let $\mathbb{S}^{n}$ be the set of $n$-by-$n$ real symmetric matrices and $\mathbb{S}^{n}_{+}\subset\mathbb{S}^{n}$ be the cone of positive semidefinite matrices. We are given a matrix $\Sigma=\left[\Sigma_{i,j}\right]_{i,j}\in\mathbb{S}^{n}$ and a number $1\leq k\leq n$. Our task is to decompose $\Sigma$ as the sum of $L=\left[L_{i,j}\right]_{i,j}\in\mathbb{S}^{n}_{+}$ of rank at most equal $k$ and a sparse matrix $S=\left[S_{i,j}\right]_{i,j}\in\mathbb{R}^{n\times n}$ in some optimal way. Observe that the matrix $S$ is also a symmetric matrix. It is well-known that the matrix $L$ can be represented as $L=MM^{T}$, where $M=\left[M_{i,j}\right]_{i,j}\in\mathbb{R}^{n\times k}$; thus $\Sigma=MM^{T}+S$.

For practical purposes, we shall represent every symmetric $n$-by-$n$ matrix by a vector of dimension $r:=n(n+1)/2$; formally, we define the linear operator

[TABLE]

The operator $h$ is obviously invertible, and its inverse shall be denoted by $h^{-1}$. Similarly, every vector of dimension $nk$ shall be represented by a $n$-by-$k$ matrix by the linear operator $g:\mathbb{R}^{nk}\to\mathbb{R}^{n\times k}$, which maps $v=(v_{1},\ldots,v_{nk})^{T}$ to $g(v):=W\equiv[w_{i,j}]_{i,j}$ with $w_{i,j}=v_{(i-1)k+j}$, in a kind of row-after-row filling of $W$ by $v$. This operator has clearly an inverse $g^{-1}$.

We construct a neural network with $n(n+1)/2$ inputs and $nk$ outputs; these outputs are meant to represent the coefficients of the matrix $M$ with whom we shall construct the rank $k$ matrix $L$ in the decomposition of the input matrix $\Sigma$. However, we do not use this neural network in its feed-forward mode as a heuristic to compute $M$ from an input $\Sigma$; we merely use the neural network framework as a way to parametrize a tentative solution $M$ to our decomposition problem.

We construct our neural network with $m$ layers of respectively $\ell_{1},\ldots,\ell_{m}$ neurons, each with the same activation function $\sigma:\mathbb{R}\to[-1,1]$. We assume that the first and the second derivative of $\sigma$ are uniformly bounded from above by the constants $\sigma^{\prime}_{\textrm{max}}>0$ and $\sigma^{\prime\prime}_{\textrm{max}}$, respectively. We let $\ell_{0}:=n(n+1)/2$ and $\ell_{m+1}:=nk$. In accordance with the standard architecture of fully connected multi-layered neural networks, for $0\leq u\leq m$ we let $A^{(u)}=\left[A^{(u)}_{i,j}\right]_{i,j}\in\mathbb{R}^{\ell_{u+1}\times\ell_{u}}$ be layer $u$’s weights,

$b^{(u)}=\left[b^{(u)}_{i}\right]_{i}\in\mathbb{R}^{\ell_{u+1}}$ be its bias, and for all $v\in\mathbb{R}^{\ell_{u}}$

[TABLE]

For each $1\leq i\leq m$ and each $v\in\mathbb{R}^{\ell_{i}}$, we write $\sigma^{(i)}(v):=\left(\sigma(v_{1}),\ldots,\sigma(v_{\ell_{i}})\right)^{T}$. We denote the parameters $\Theta:=(A^{(0)},b^{(0)},\dots,A^{(m)},b^{(m)})$ and define the $m$-layered neural network by the function $\mathcal{N}^{\Theta}_{m}$ from $\mathbb{R}^{\ell_{0}}$ to $\mathbb{R}^{\ell_{m+1}}$ for which $\mathcal{N}^{\Theta}_{m}(x)$ equals

[TABLE]

We therefore have to specify $\mathbb{L}_{m}:=\sum_{u=0}^{m}\ell_{u}\ell_{u+1}+\ell_{u+1}$ many parameters to describe the neural network $\mathcal{N}^{\Theta}_{m}$ completely.

Now, we are ready to define the cost function to minimize. We write the $1$-norm of some $X\in\mathbb{R}^{n\times n}$ as $||X||_{1}:=\sum_{i,j=1}^{n}|X_{i,j}|$. Our objective function is, for a given $\Sigma\in\mathbb{S}^{n}$, the function

[TABLE]

Since M=g\big{(}\mathcal{N}^{\Theta}_{m}\left(h(\Sigma)\right)\!\big{)} is our tentative solution to the matrix decomposition problem, this objective function consists in minimizing $||MM^{T}-\Sigma||_{1}=||S||_{1}$ over the parameters $\Theta$ that define $M$.

As the function $\Theta\mapsto\varphi_{\mathrm{obj}}(\Theta)$ is neither differentiable nor convex, we do not have access neither to its gradient nor to a subgradient. We shall thus approximate it by

[TABLE]

where $\mu:\mathbb{R}\to[0,\infty)$ is a smooth approximation of the absolute value function with a derivative uniformly bounded by $1$ and its second derivative bounded by $\mu^{\prime\prime}_{\max}$. A widely used example of such a function is given by

[TABLE]

where $\varepsilon$ is a small positive constant. With this choice for $\mu$, we have $\mu^{\prime\prime}_{\max}=1/\varepsilon$. Another example, coming from the theory of smoothing techniques in convex optimization [17], is given by $\mu(t):=\varepsilon\ln(2\cosh(t/\varepsilon)),$ also with $\mu^{\prime\prime}_{\max}=1/\varepsilon$.

We apply a gradient method to minimize the objective function ${\varphi}$, whose general scheme can be written as follows.

[TABLE]

The norm we shall use in the sequel is a natural extension of the standard Frobenius norm to finite lists of matrices of diverse sizes (the Frobenius norm of a vector coinciding with its Euclidean norm). Specifically, for any $\gamma\in\mathbb{N}_{0}$, $m_{1},\ldots,m_{\gamma},n_{1},\ldots,n_{\gamma}\in\mathbb{N}_{0}$, and $(X^{1},\ldots,X^{\gamma})\in\mathbb{R}^{m_{1}\times n_{1}}\times\dots\times\mathbb{R}^{m_{\gamma}\times n_{\gamma}}$, we let

[TABLE]

This norm is merely the standard Euclidean norm of the vector obtained by concatenating all the columns of $X^{1},\ldots,X^{\gamma}$.

Since the objective function in (9) is non-convex, this method can only realistically converge to one of its stationary point or to stop close enough from one, that is, at point $\Theta^{*}$ for which $||\nabla\varphi(\Theta^{*})||$ is smaller than a given tolerance. The complexity of many variants of this method can be established if the function ${\varphi}$ has a Lipschitz continuous gradient (see [8] and references therein). We have the following convergence result.

Theorem II.1**.**

Let $\Sigma\in\mathbb{S}^{n}$ and assume that there exists $D>0$ such that the sequence $(\Theta_{j})_{j\in\mathbb{N}_{0}}$ of parameters constructed in (10) satisfies

[TABLE]

Then, the gradient of the function $\varphi$ defined in (9) is Lipschitz continuous on $\mathcal{D}:=\{\Theta\in\mathbb{R}^{\mathbb{L}_{m}}\colon\|\Theta\|\leq D\}$ with Lipschitz constant bound $L_{m}$ satisfying for $m\geq 1$:

[TABLE]

where $\mathcal{C}_{m}$ is a constant that only depends polynomially on $\sigma^{\prime}_{\max},\sigma^{\prime\prime}_{\max},\mu^{\prime\prime}_{\max}$, with powers in $O(m)$ and $\ell_{\max}:=\max\{\ell_{1},\ldots,\ell_{m}\}$. When $m=1$, we have an alternative Lipschitz constant bound, more favourable for large $D$:

[TABLE]

As a consequence, if for the gradient method (10) there exists a constant $K>0$ such that for all $j\geq 0$

[TABLE]

then for every $N\in\mathbb{N}$ we have that

[TABLE]

where $\varphi^{*}:=\min_{\Theta\in\mathcal{D}}\varphi(\Theta)$. In particular, for every tolerance level $\varepsilon>0$ we have

[TABLE]

Remark II.2** (Choosing $(h_{j})$).**

Notice that the condition (13) in Theorem II.1 imposed on the gradient method, or more precisely on the step-size strategy $(h_{j})$, is not very restrictive. We provide several examples that are frequently used.

*The sequence $(h_{j})$ is chosen in advance, independently of the minimization problem. This includes, e.g., the common constant step-size strategy $h_{j}=h$ or $h_{j}=\tfrac{h}{\sqrt{j+1}}$ for some constant $h>0$. With these choices, one can show that (13) is satisfied for $K=1$. * 2. 2.

The Goldstein-Armijo rule, in which, given $0<\alpha<\beta<1$, one needs to find $(h_{j})$ such that

[TABLE]

This strategy satisfies (13) with $K=2\alpha(1-\beta)$.

We refer to [18, Section 1.2.3] and to [19, Chapter 3] for further details and more examples.

Remark II.3** (On Assumption (12)).**

The convergence rate (14) obtained in Theorem II.1 relies fundamentally on the Lipschitz property of the gradient of the (approximated) objective function $\varphi$ of the algorithm in (10). However, due to its structure, we see that the global Lipschitz property of $\nabla\varphi$ fails already for a single-layered neural network, as it grows polynomially of degree 4 in the parameters; see also Section IV. Yet, it is enough to ensure the Lipschitz property of $\nabla\varphi$ on the domain of the sequence of parameters $(\Theta_{j})_{j\in\mathbb{N}_{0}}$ generated by the algorithm in (10), which explains the significance of assumption (12). Nevertheless, assumption (12) is not very restrictive as one might expect that the algorithm (10) converges and hence automatically forces assumption (12) to hold true. Empirically, we verify in Subsection III-D that this assumption holds for our algorithm when used with our two non-synthetic data sets.

Remark II.4** (Polynomiality of our method).**

While the second part of Theorem II.1 is standard in optimization (see, e.g., in [18, Section 1.2.3]), we notice that for a fixed depth $m$ of the neural network the constant $L$ in the rate of convergence of the sequence $(\min_{1\leq j\leq N}\|\nabla\varphi(\Theta_{j})\|)$ grows polynomially in the parameters $\ell_{\max}:=\max\{\ell_{1},\dots,\ell_{m}\}$, $n$, and $k$. These parameters describe the dimensions of the input, the output, and the hidden layers of the neural network. A rough estimate shows that

[TABLE]

Remark II.5** (A simplified version of (14)).**

Note that, since the function $\varphi$ is nonnegative, we have $\varphi^{*}\geq 0$, so that (14) can be simplified by

[TABLE]

Remark II.6** (Choice of activation function).**

We require the activation function $\sigma\colon\mathbb{R}\to[-1,1]$ to be a non-constant, bounded, and smooth function. The following table provides the most common activation functions satisfying the above conditions.

We provide the proof of Theorem II.1 in Section IV.

III Numerical results

III-A Numerical results based on simulated data

We start our numerical tests111We refer to https://github.com/CalypsoHerrera/lrs-cov-nn for the codes of our simulations. with a series of experiments on artificially generated data. We construct a collection of $n$-by-$n$ positive semidefinite matrices $\Sigma$ that can be written as $\Sigma=L_{0}+S_{0}$ for a known matrix $L_{0}$ of rank $k_{0}\leq n$ and a known matrix $S_{0}$ of given sparsity $s_{0}$. We understand by sparsity the number of null elements of $S_{0}$ divided by the number of coefficients of $S_{0}$; when a sparse matrix is determined by an algorithm, we consider that every component smaller in absolute value than $\epsilon=0.01$ is null. To construct one matrix $L_{0}$, we first sample $nk_{0}$ independent standard normal random variables that we arrange into an $n$-by-$k_{0}$ matrix $M$. Then $L_{0}$ is simply taken as $MM^{T}$. To construct a symmetric positive semidefinite sparse matrix $S_{0}$ we first select uniformly randomly a pair $(i,j)$ with $1\leq i<j\leq n$. We then construct an $n$-by-$n$ matrix $A$ that has only four non-zero coefficients: its off-diagonal elements $(i,j)$ and $(j,i)$ are set to a number $b$ drawn uniformly randomly in $[-1,1]$, whereas the diagonal elements $(i,i)$ and $(j,j)$ are set to a number $a$ drawn uniformly randomly in $[|b|,1]$. This way, the matrix $A$ is positive semidefinite. The matrix $S_{0}$ is obtained by summing such matrices $A$, each corresponding to a different pair $(i,j)$, until the desired sparsity $s_{0}$ is reached.

Given an artificially generated matrix $\Sigma=L_{0}+S_{0}$, where $L_{0}$ has a prescribed rank $k_{0}$ and $S_{0}$ a sparsity $s_{0}$, we run our algorithm to construct a matrix $M\in\mathbb{R}^{n\times k}$. With $L:=MM^{T}$ and $S:=\Sigma-L$, we determine the approximated rank $r(L)$ of $L$ by counting the number of eigenvalues of $L$ that are larger than $\epsilon=0.01$. We also determine the sparsity $s(S)$ as specified above, by taking as null every coefficient smaller than $\epsilon=0.01$ in absolute value. We compute the discrepancy between the calculated low-rank part $L$ and the correct one $L_{0}$ by $\textrm{rel.error}(L):=||L-L_{0}||_{F}/||L_{0}||_{F}$ and between $S$ and the true $S_{0}$ by $\textrm{rel.error}(S):=||S-S_{0}||_{F}/||S_{0}||_{F}$. Table II reports the average of these quantities over ten runs of our algorithm DNN (short for Deep Neural Network), as well as their standard deviation (in parenthesis). We carried our experiments on various values for the dimension $n$ of the matrix $\Sigma$, for the given rank $k_{0}$ of $L_{0}$, for the given sparsity $s_{0}$ of $S_{0}$ and for the chosen forced (upper bound for the) rank $k$ in the construction of $L$ introduced in Section II.

When choosing $n=100$, our algorithm unsurprisingly achieves the maximal rank $k$ for the output matrix $L$, unlike IALM, Non-Convex RPCA, and FPCP. The sparsities are comparable when $n=100$, even though the different methods have different strategies to sparsify their matrices; some methods, such as FPCP, apply a shrinkage after optimization by replacing every matrix entries $S_{ij}$ by $\operatorname{sign}(S_{ij})[|S_{ij}|-1/\sqrt{n}]^{+}$. By forcing sparsity, FPCP makes its output violate Equation (7), so the sparsity of its output and ours might not be comparable. We however display $s(S)$ and its relative error $\textrm{rel.error}(S)$ even for that method, to give some insight on how the algorithms behave.

When the given rank $k_{0}$ coincides the forced rank $k$, our algorithm achieves a much higher accuracy than all the other algorithms for small $n$. When $n$ gets larger, our algorithm still compares favourably with the other ones, except on a few outlier instances. Of course, when there is a discrepancy between $k_{0}$ and $k$, our algorithm cannot recover $k_{0}$. Nevertheless, the relative error $\textrm{rel.error}(L)$ compares favorably with the other methods, especially when $k_{0}>k$. We acknowledge however that in circumstances where one needs to minimize the rank of $L$, e.g., to avoid overfitting past data, the forced rank $k$ can only be considered as the maximum rank of $L$ returned by the algorithm; in such case, the PCP, IALM, and FPCP algorithms could be more appropriate.

Various network architectures and corresponding activation functions have been tested. They only marginally influence the numerical results.

III-B Application on a five hundred S&P500 stocks portfolio

In this section, we evaluate our algorithm on real market data and compare it the other algorithms we have selected to demonstrate its capability also when the low-rank plus sparse matrix decomposition is not known. A natural candidate for our experiment is the correlation matrix of stocks in the S&P500, due to its relatively large size and the abundant, easily available data. Five hundred S&P500 stocks were part of the index between 2017 and 2018. To make the representation more readable, we have sorted these stocks in eleven sectors according to the global industry classification standard222First, we have those belonging to the energy sector, second, those from the materials sector, then, in order, those from industrials, real estate, consumer discretionary, consumer staples, health care, financials, information technology, communication services, and finally utilities. We may notice that utilities seem almost uncorrelated to the other sectors, and that real estate and health care present a significantly lower level of correlation to the other sectors than the rest.. We have constructed the correlation matrix $\Sigma$ from the daily returns of these 500 stocks during 250 consecutive trading days (see Figure 1). As the data used to construct $\Sigma$ are available at an identical frequency, the matrix $\Sigma$ is indeed positive semidefinite, with 146 eigenvalues larger than $10^{-10}$. The 70 largest eigenvalues account for 90% of $\Sigma$’s trace, that is, the sum of all its 500 eigenvalues.

In Figure 1, we display the resulting matrices $L$ and $S$ for all algorithms with respect to the same input $\Sigma$. In our method and in RPCA-GD, we have set the rank of $L$ to $k=3$. Coincidentally, we have also obtained a rank of $3$ with FPCP. Among the three output matrices $L$, the one returned by our method matches the input more closely as it contains more relevant eigenspaces. The other algorithms have transferred this information to the sparse matrix $S$. Note that the scale of values for the correlation matrix ranges between $-0.2$ and $0.9$. The ranks of $L$ obtained with the two other algorithms are $61$ for IALM and $5$ for Non-Convex RPCA, showing the difficulty of tuning these methods to obtain some desired rank. The matrix $S$ from DNN would normally be slightly less sparse than the matrix $S$ of FPCP, as FPCP applies shrinkage. However, for visualization and comparison purposes, shrinkage is also applied on $S$ returned by DNN algorithm in Figure 1.

III-C Application on real estate return

We have computed the low-rank plus sparse decomposition of the real estate return matrix for $44$ countries333The countries are ordered by continent and subcontinents: Western Europe (Belgium, Luxembourg, Netherlands, France, Germany, Switzerland, Austria, Denmark, Norway, Sweden, Finland, United Kingdom, Ireland, Italy, Spain, Portugal), Eastern Europe (Croatia, Estonia, Latvia, Lithuania, Russia, Poland, Bulgaria, Hungary, Romania, Slovak Republic, Czech Republic), Near East (Turkey, Saudi Arabia), Southern America (Brazil, Chile, Colombia, Peru, Mexico), Northern America (United States, Canada), Eastern Asia (India, China, Hong Kong, Singapore, Japan) Oceania (Australia, New Zealand) and Africa (South Africa).. The correlation matrix contains $88$ returns, alternating the residential returns and the corporate returns of each country444 We thank Eric Schaanning for providing us this correlation matrix.; see Figure 2. We impose the rank of the output matrix to be equal to $3$. Similar to the previous section, the correlation color scale in Figure 2 is cropped between $-0.5$ and $0.5$ for a better visualization. The sparse matrix of FPCP is normally sparser than the one returned by our DNN algorithm. However in Figure 2, and for a fair comparison, shrinkage is also applied to the matrix $S$ returned by our algorithm. The low-rank matrices $L$ exhibit a variety of ranks from low (1 for NC-RPCA) to high (14 for IALM), indicating the difficulty of tuning the hyperparameters of these methods to a desired rank.

In Figure 3, we plot in the first line the eigenvalues of the matrix $L$ returned by FPCP and DNN, as well as the eigenvalues of the original matrix $\Sigma$, and for all the algorithms simultaneously in the second line. In the left figure, we display the first $17$ eigenvalues of the matrix $L$ where the forced rank is set to $k=15$. In the right figure, we plot the first $50$ eigenvalues where the forced rank is set to $k=88$. Notice that the input matrix $\Sigma$ has some negative eigenvalues. This phenomenon can happen in empirical correlation matrices when the data of the different variables are either not sampled over the same time frame or not with the same frequency; we refer to [11] for a further discussion on this issue. Our DNN algorithm and RPCA-GD, by setting $L:=MM^{T}$, avoid negative eigenvalues, although the original matrix $\Sigma$ is not positive semidefinite. In contrast, the other algorithms might output a non-positive semidefinite matrix.

III-D Empirical verification of bounded parameters

To verify empirically our assumption (12) in Theorem II.1 that the parameters $(\Theta_{j})_{j\in\mathbb{N}_{0}}$ generated by our algorithm (10) remain in a compact set, we plotted in Figure 4 the running maximum $\max_{0\leq j\leq J}\|\Theta_{j}\|$ as a function of the number of iterations $J$ for both examples on the S&P500 and the real estate data used in the previous section. For both cases, we observe, as desired, that the running maximum $\max_{0\leq j\leq J}\|\Theta_{j}\|$ converges, which means that at least empirically, $(\Theta_{j})_{j\in\mathbb{N}_{0}}$ remains in a compact set.

IV Proof of convergence

A vast majority of first-order methods for minimizing locally a non-convex function with provable convergence rate are meant to minimize $L$-smooth functions, that is, differentiable functions with a Lipschitz continuous gradient. Also, the value of the Lipschitz constant with respect to a suitable norm plays a prominent role in this convergence rate (see [8] and references therein). As a critical step in the convergence proof for the minimization procedure of the function ${\varphi}$, we compute carefully a bound on the Lipschitz constant of its gradient, also called the smoothness constant below.

For establishing this bound, we use in a critical way the recursive compositional nature of the neural network that represents $M$: we first consider the case of a neural network with a single layer. We compute explicitly the derivative of the corresponding objective function $\varphi(\Theta)$ with respect to each parameter in Lemmas IV.1 and IV.2. Then, using an elementary lemma on the smoothness constant of a composition and a product of functions, we carefully derive an upper bound on the derivative of a single-layered network output with respect to the bias on the input in Lemma IV.4. This result represents the key element of our proof. Lemma IV.5 completes the picture by computing such a bound for derivatives with respect to other coefficients. Lemma IV.6 merges these results and provides a bound on the overall smoothness constant for the case where $M$ is parametrized by a single-layer neural network. Finally, we deduce inductively a bound on this Lipschitz constant for a multi-layered neural network in the proof of Theorem II.1.

So, we consider in the beginning of this section a neural network with a single layer of $\ell$ neurons. These neurons have each $\sigma:\mathbb{R}\to[-1,1]$ as activation function, and $\sigma$ has its first and the second derivative uniformly bounded by $\sigma^{\prime}_{\textrm{max}}>0$ and $\sigma^{\prime\prime}_{\textrm{max}}$, respectively. We let $r:=n(n+1)/2$, $A=\left[A_{i,j}\right]_{i,j}\in\mathbb{R}^{\ell\times r}$ be the weights, and $b=\left[b_{i}\right]_{i}\in\mathbb{R}^{\ell}$ be the bias on the input, and let

[TABLE]

The coefficients on the output are denoted by $C=\left[C_{i,j}\right]_{i,j}\in\mathbb{R}^{nk\times\ell}$ and the bias by $d=\left[d_{i}\right]_{i}\in\mathbb{R}^{nk}$. As above, we define

[TABLE]

With $\tilde{\sigma}:\mathbb{R}^{\ell}\to\mathbb{R}^{\ell}$, $u\mapsto\tilde{\sigma}(u)=\left(\sigma(u_{1})\cdots\sigma(u_{\ell})\right)^{T}$, and $\Theta:=(A,b,C,d)$ the single layer neural network $\mathcal{N}^{\Theta}:\mathbb{R}^{r}\to\mathbb{R}^{nk}$ is then the composition of these three functions, that is: $\mathcal{N}^{\Theta}=f_{2}^{C,d}\circ\tilde{\sigma}\circ f_{1}^{A,b}$. For a given $\Sigma\in\mathbb{S}^{n}$, our approximated objective function with respect to the above single-layer neural network is defined by

[TABLE]

where $\mu:\mathbb{R}\to\mathbb{R}$ is a smooth approximation of the absolute value function with a derivative uniformly bounded by $1$ and its second derivative bounded by $\mu^{\prime\prime}_{\max}$.

As announced above, we start our proof by computing the partial derivatives of $\tilde{\varphi}$. For abbreviating some lengthy expressions, we use the following shorthand notation throughout this section. We fix $\Sigma\in\mathbb{S}^{n}$ and let for $1\leq i,j\leq n$,

[TABLE]

Lemma IV.1**.**

Let $\tilde{\varphi}$ be the function defined in (15) and $1\leq\iota\leq\ell$. Then

[TABLE]

Moreover, for every $1\leq\iota\leq\ell$, $1\leq\eta\leq r$ we have

[TABLE]

Proof.

Let $1\leq\iota\leq\ell$. By definition of $\tilde{\varphi}$, we have

[TABLE]

As $\mathcal{N}^{\Theta}\equiv f_{2}^{C,d}\circ\tilde{\sigma}\circ f_{1}^{A,b}$, we have

[TABLE]

Observe that

[TABLE]

The third order tensor $\nabla_{X}g(X)$, as an $nk$-dimensional vector of $n$-by-$k$ matrices, has for $((i-1)k+j)$-th element the matrix whose only nonzero element is a $1$ at position $(i,j)$, namely

[TABLE]

Plugging these in (18), we get

[TABLE]

so that $\frac{\partial}{\partial b_{\iota}}\tilde{\varphi}(\Theta)$ equals

[TABLE]

and the first part is proved. For the second part, note that for every $1\leq\eta\leq r$ we have

[TABLE]

which coincides with $[h(\Sigma)]_{\eta}\frac{\partial f_{1}^{A,b}(h(\Sigma))}{\partial b_{\iota}}$ by (22). Hence, using the same derivations as in (17) and (18), we get

[TABLE]

∎

Lemma IV.2**.**

Let $\tilde{\varphi}$ be the function defined in (15), fix $1\leq\alpha\leq n$ and $1\leq\beta\leq k$, and set $\nu:=(\alpha-1)k+\beta$. Then

[TABLE]

Moreover, for every $1\leq\iota\leq\ell$, we have

[TABLE]

Proof.

Note that $\frac{\partial\tilde{\varphi}(\Theta)}{\partial d_{\nu}}$ equals

[TABLE]

As $\nu=(\alpha-1)k+\beta$, the same calculation as in (22) and (23) ensures that

[TABLE]

which is a $n$-by-$k$ matrix. Hence, $\frac{\partial\tilde{\varphi}(\Theta)}{\partial d_{\nu}}$ equals

[TABLE]

Since $\omega_{i,j}=\omega_{j,i}$ for every $1\leq i,j\leq n$ we finally get

[TABLE]

which proves the first part. For the second part, we fix $1\leq\iota\leq\ell$. In the same manner as we established (24), we can use (23) to see that

[TABLE]

and thereby, with (26), to conclude that

[TABLE]

∎

One of the key tools in the derivation of bounds on the Lipschitz constant of $\tilde{\varphi}$ is the following elementary lemma. It shows how to infer the Lipschitz constant of some functions from the Lipschitz constant of other, simpler Lipschitz-continuous functions.

Lemma IV.3**.**

Let $\phi_{1}:\mathbb{R}^{m}\to\mathbb{R}$ be Lipschitz continuous with Lipschitz constant $L_{1}$. Let $\phi_{2}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ be a function for which there exists a some $L_{2}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to(0,\infty)$ such that for all $x,y\in\mathbb{R}^{n}$ we have

[TABLE]

Assume that $\phi_{1}\circ\phi_{2}$ is bounded by a constant $B_{12}$. Finally, let $\phi_{3}:\mathbb{R}^{n}\to\mathbb{R}$ be a function for which

$|\phi_{3}(y)|\leq B_{3}(y)$* for all $y\in\mathbb{R}^{n}$ for some positive function $B_{3}\colon\mathbb{R}^{n}\to(0,\infty)$;* 2. 2.

there exist three functions $L_{31}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to(0,\infty)$, $L_{32}:\mathbb{R}^{n}\times\mathbb{R}^{n}\to(0,\infty)$, and $f:\mathbb{R}^{n}\to\mathbb{R}^{p}$ such that for all $x,y\in\mathbb{R}^{n}$ we have

[TABLE]

Then, for every $x,y\in\mathbb{R}^{n}$, the function $\Phi:=(\phi_{1}\circ\phi_{2})\phi_{3}$ verifies

[TABLE]

Proof.

For all $x,y\in\mathbb{R}^{n}$, we can write

[TABLE]

It remains to use (28) on the last term to conclude. ∎

As defined in (11), the norm we use for a finite list of matrix of different sizes is the Euclidean norm of the vector constituted by their Frobenius norms. For symmetric matrices $X$, we also use a dedicated norm defined as $||X||_{S}:=||h(X)||$, where $h:\mathbb{S}^{n}\to\mathbb{R}^{n(n+1)/2}$ is the function defined in (8). Note also that $2||X||_{S}^{2}=||X||_{F}^{2}+||\textrm{diag}(X)||_{F}^{2}$.

The next lemma is a key step in determining a bound on the smoothness constant of our objective function $\tilde{\varphi}$.

Lemma IV.4**.**

Let $\tilde{\varphi}$ be the function defined in (15) and let $D>0$ be any constant. We define

[TABLE]

The function

[TABLE]

is Lipschitz continuous on $\mathcal{D}$ with a constant $L_{b}$ for which

[TABLE]

where $L_{Z}:=\sqrt{1+\|\Sigma\|_{S}^{2}}$ and $\mathcal{C}_{b}$ is a constant that only depends polynomially on $\sigma^{\prime}_{\max},\sigma^{\prime\prime}_{\max},\mu^{\prime\prime}_{\max}$.

Proof.

We divide the proof into several steps, each of which establishing that some function is bounded and/or Lipschitz continuous. We set $\Theta=(A,b,C,d)$ and $\bar{\Theta}=(\bar{A},\bar{b},\bar{C},\bar{d})$ in $\mathcal{D}$. We make use abundantly of the shorthand notation defined in (16), adopting the notation $\bar{\omega}_{ij}$, $\bar{X}$, etc…when $\Theta$ is replaced by $\bar{\Theta}$.

Step 1. Let $1\leq\iota\leq\ell$. In this first step, we focus on $Z_{\iota}=\left[f_{1}^{A,b}(h(\Sigma))\right]_{\iota}$. The Cauchy-Schwarz Inequality ensures that for every $(A,b),(\bar{A},\bar{b})\in\mathbb{R}^{\ell\times r}\times\mathbb{R}^{\ell}$ we have

[TABLE]

Step 2. Let us apply Lemma IV.3 with $\phi_{1}\equiv\sigma$, $\phi_{2}\equiv\left[f_{1}^{\cdot,\cdot}(h(\Sigma))\right]_{\iota}$, and $\phi_{3}\equiv 1$. We have immediately $L_{31}=L_{32}=0$, $B_{3}=1$, $L_{1}=\sigma^{\prime}_{\max}$, and, by Step $1$, $L_{2}=L_{Z}$. Therefore,

[TABLE]

Note also that $\|\tilde{\sigma}(\bar{Z})\|\leq\sqrt{\ell}$ because $|\sigma\left(\bar{Z}_{\iota}\right)|\leq 1$ by assumption on $\sigma$.

Step 3. Similarly, we get

[TABLE]

Also, $\sigma^{\prime}\left(Z_{\iota}\right)\leq\sigma^{\prime}_{\max}$ by assumption.

Step 4. Let us fix $1\leq\nu\leq nk$. We focus in this step on $X_{\nu}=\left[f_{2}^{{C,}{d}}\circ\tilde{\sigma}\circ f_{1}^{{A},{b}}(h(\Sigma))\right]_{\nu}=\left[f_{2}^{{C,}{d}}(Y)\right]_{\nu}$, with $Y=\tilde{\sigma}(Z)=\tilde{\sigma}\circ f_{1}^{{A},{b}}(h(\Sigma))$. Observe first that

[TABLE]

Squaring both sides of this inequality and using Cauchy-Schwarz, we get

[TABLE]

with

[TABLE]

This last inequality is also ensured by Step $2$, where we showed that $||\bar{Y}||\leq\sqrt{\ell}$.

In addition to providing an estimate for its Lipschitz constant, we can also obtain a bound for this function. From Step $2$, we get

[TABLE]

We obtain with this step the Lipschitz constant of one output of a single-layer neural network. Would our objective function be the standard least square loss function to minimize when training a neural network, the remaining of our task would have been vastly simpler. The specific intricacies of our objective function (15) ask for a more involved analysis.

Step 5. Let $1\leq i,j\leq n$. In this step, we deal with the function

[TABLE]

We have

[TABLE]

Let us bound these four sums using (29) and (31). First,

[TABLE]

The sum on the right-hand side can be conveniently rewritten using the vector-to-matrix operator $g$.

[TABLE]

Similarly, we can write $\sum_{s=1}^{k}\|d_{(i-1)k+s}\|^{2}=\|g(d)_{i,:}\|^{2}$, so that

[TABLE]

Second, with $\Theta_{I(\nu)}:=(A,b,C_{\nu,:},d_{\nu})$ and writing $\nu(s):=(j-1)k+s$ for $1\leq s\leq k$, we have

[TABLE]

We can use the crude bound

[TABLE]

to get

[TABLE]

by (33). We can plug our two estimates into our bound for $|\omega_{i,j}-\bar{\omega}_{i,j}|^{2}$ to get

[TABLE]

with

[TABLE]

where, $\kappa:=(\ell+1)/\left(\sigma^{\prime}_{\max}L_{Z}\right)^{2}$.

Step 6. In this step, we fix $1\leq i,j\leq n$ and $1\leq\iota\leq\ell$ and focus on $\mu^{\prime}\left({\omega}_{i,j}\right)\sigma^{\prime}({Z}_{\iota})$. We therefore seek to apply Lemma IV.3 with $\phi_{1}\equiv\mu^{\prime}$, $\phi_{2}\equiv{\omega}_{i,j}$, and $\phi_{3}\equiv\sigma^{\prime}({Y}_{\iota})$. The assumptions on $\mu$ allow us to take $L_{1}=\mu^{\prime\prime}_{\max}$ and $B_{12}=1$. Step $5$ shows that $L_{2}$ can be taken as $L_{\omega,i,j}$. Step $3$ allows us to pick $L_{31}=\sigma^{\prime\prime}_{\max}L_{Z}$ with $L_{32}=0$ and $B_{3}=\sigma^{\prime}_{\max}$ by our assumptions on $\sigma$.

We obtain that

[TABLE]

with

[TABLE]

Evidently, we also have $\mu^{\prime}\left({\omega}_{i,j}\right)\sigma^{\prime}({Z}_{\iota})\leq\sigma^{\prime}_{\max}$.

Step 7. The formula for $\frac{\partial\tilde{\varphi}(\Theta)}{\partial b_{\iota}}$ obtained in Lemma IV.1 appears as a weighted sum of sums. Step $6$ did focus on the weights of this combination. In this step, we deal with these sums. We fix again $1\leq i,j\leq n$ and $1\leq\iota\leq\ell$ and define

[TABLE]

We define similarly $\bar{V}_{i,j,\iota}$ according to our usual convention to replace $\Theta$ by $\bar{\Theta}$. Using the same argument as in Step $5$, we can write

[TABLE]

We can bound the sums involving $X$ or $\bar{X}$ as in Step $5$. The sums with $C,\bar{C}$ can we rewritten using (33); in particular,

[TABLE]

We obtain

[TABLE]

with

[TABLE]

We can also obtain an upper bound to \big{|}V_{i,j,\iota}\big{|} due to

[TABLE]

by using Cauchy-Schwarz, (33), and a bound obtained in Step $5$.

Step 8. We start this final step by applying Lemma IV.3 to the function $\Upsilon_{i,j,\iota}:=\mu^{\prime}(\omega_{i,j})\sigma^{\prime}(Z_{\iota})\left(V_{i,j,\iota}+V_{j,i,\iota}\right)$ by taking $\phi_{1}\equiv id_{\mathbb{R}}$, $\phi_{2}\equiv\mu^{\prime}(\omega_{i,j})\sigma^{\prime}(Z_{\iota})$, and $\phi_{3}\equiv V_{i,j,\iota}+V_{j,i,\iota}$. Evidently $L_{1}=1$. In Step $6$, we have shown that we can take $B_{12}=\sigma^{\prime}_{\max}$ and $L_{2}=L_{\mu^{\prime}\sigma^{\prime},i,j}$. By Step $7$, we can set $B_{3}=B_{CX,i,j,\iota}+B_{CX,j,i,\iota}$, $L_{31}=L_{CdX,i,j,\iota}+L_{CdX,j,i,\iota}$, and $L_{32}=0$. We obtain that

[TABLE]

Since in view of Lemma IV.4, $\frac{\partial\tilde{\varphi}(\Theta)}{\partial b_{\iota}}$ is the sum of $\Upsilon_{i,j,\iota}$ over $i,j$, we obtain immediately the Lipschitz continuity of this partial derivative. Squaring the above bound and using Cauchy-Schwarz on the 6-terms right-hand side, we get

[TABLE]

It remains to estimate these six sums. For the first one (and the second one by symmetry), we can proceed as follows.

[TABLE]

To estimate the third and fourth sum, notice that

[TABLE]

For the last two terms, we can write

[TABLE]

Therefore, using all the estimates for the six terms, we arrive at

[TABLE]

We conclude that

[TABLE]

where $\mathcal{C}_{b}\equiv\mathcal{C}_{b}(\sigma^{\prime}_{\max},\sigma^{\prime\prime}_{\max},\mu^{\prime\prime}_{\max})$ is a constant that only depends polynomially on $\sigma^{\prime}_{\max},\sigma^{\prime\prime}_{\max},\mu^{\prime\prime}_{\max}$. ∎

In the next lemma, we determine a bound for the Lipschitz constant of $\frac{\partial\tilde{\varphi}(\Theta)}{\partial d_{\nu}}$. Lemmas IV.1 and IV.2 show how to deduce from the previous lemma a bound on the Lipschitz constant of $\frac{\partial\tilde{\varphi}(\Theta)}{\partial A_{\iota,\eta}}$ and of $\frac{\partial\tilde{\varphi}(\Theta)}{\partial C_{\nu,\iota}}$ from the next lemma. We use the same set $\mathcal{D}$ as in Lemma IV.4.

Lemma IV.5**.**

Let $\tilde{\varphi}$ be the function defined in (15) and let us consider the set $\mathcal{D}$ defined in the statement of Lemma IV.4 for some $D>0$. The function

[TABLE]

is Lipschitz continuous on $\mathcal{D}$ with a constant $L_{d}$ satisfying

[TABLE]

where $L_{Z}:=\sqrt{1+\|\Sigma\|_{S}^{2}}$ and $\mathcal{C}_{d}$ is a constant that depends polynomially on $\sigma^{\prime}_{\max}$ and $\mu^{\prime\prime}_{\max}$.

Proof.

We set $\Theta=(A,b,C,d)$ and $\bar{\Theta}=(\bar{A},\bar{b},\bar{C},\bar{d})$ two distinct points in $\mathcal{D}$. Let us fix $1\leq\nu\leq nk$ and set $\alpha$ and $\beta$ to be the unique numbers for which $\nu=(\alpha-1)k+\beta$, $1\leq\alpha\leq n$ and $1\leq\beta\leq k$. We computed in Lemma IV.2

[TABLE]

with $\omega_{\alpha,j}$ and $X$ as defined in (16).

Let $1\leq j\leq n$ and set $\eta(j):=(j-1)k+\beta$. To obtain a bound on the Lipschitz constant of $\Gamma_{\alpha,\beta,j}:=\mu^{\prime}(\omega_{\alpha,j})X_{\eta(j)}$, we apply Lemma IV.3 with $\phi_{1}\equiv\mu^{\prime}$, $\phi_{2}\equiv{\omega}_{\alpha,j}$, and $\phi_{3}\equiv X_{\eta(j)}$. By assumption, we know that $L_{1}:=\mu^{\prime\prime}_{\max}$ and $B_{12}:=\mu^{\prime}_{\max}\leq 1$. By Step $5$ in the proof of Lemma IV.4, we can take $L_{2}:=L_{\omega,\alpha,j}$, while Step $4$ ensures we can set $L_{31}:=L_{\eta(j)}$ as in (30) (note that $L_{2}$ and $L_{31}$ are functions of $\Theta$ and $\bar{\Theta}$). Finally, we can let $B_{3}:=\sqrt{\ell+1}\cdot\sqrt{||{C}_{\eta(j),:}||^{2}+|{d}_{\eta(j)}|^{2}}$ and $L_{32}=0$. We deduce that

[TABLE]

with

[TABLE]

Since

[TABLE]

the function $\frac{\partial\tilde{\varphi}(\Theta)}{\partial d_{\nu}}$ is Lipschitz continuous. Now, the Lipschitz constant of $\frac{\partial\tilde{\varphi}(\Theta)}{\partial d}$ can be estimated by

[TABLE]

For the first sum above, we can write first

[TABLE]

where this last equality, similarly to (33), follows from

[TABLE]

Therefore,

[TABLE]

For the second sum above, we have from (34):

[TABLE]

with $\kappa:=(\ell+1)/\left(\sigma^{\prime}_{\max}L_{Z}\right)^{2}$. Observe that

[TABLE]

by the classical inequality $\sum_{i}\lambda_{i}^{2}\leq\left(\sum_{i}|\lambda_{i}|\right)^{2}$. Thus,

[TABLE]

Putting everything together, we arrive at

[TABLE]

Therefore,

[TABLE]

where $\mathcal{C}_{d}$ is a constant that only depends polynomially on $\sigma^{\prime}_{\max},\mu^{\prime\prime}_{\max}$. ∎

We can now merge all our previous results to determine a bound on the Lipschitz constant of $\nabla\tilde{\varphi}$. Again, we use the set $\mathcal{D}$ as in Lemma IV.4 and $L_{Z}:=\sqrt{1+||\Sigma||_{S}^{2}}$.

Lemma IV.6**.**

The function

[TABLE]

is Lipschitz continuous with Lipschitz constant $L_{\tilde{\varphi}}$ satisfying

[TABLE]

where $\mathcal{C}_{\tilde{\varphi}}$ is a constant that only depends polynomially on $\sigma^{\prime}_{\max}$, $\sigma^{\prime\prime}_{\max}$, and $\mu^{\prime\prime}_{\max}$.

Proof.

Let $\Theta=(A,b,C,d)$ and $\bar{\Theta}=(\bar{A},\bar{b},\bar{C},\bar{d})$ two distinct points in $\mathcal{D}$. In view of Lemmas IV.1 and IV.2

[TABLE]

because $Y_{\iota}=\sigma([Ah(\Sigma)+b]_{\iota})\leq 1$ for all $1\leq\iota\leq\ell$. It remains now to use our estimates for $L_{b}$ and $L_{d}$ obtained in Lemma IV.4 and Lemma IV.5. ∎

Now that the single-layered neural network parametrization has been completely treated, we can turn our attention to a multi-layered neural network configuration. The set $\mathcal{D}$ is now the one defined in the statement of Theorem II.1, and the notation follows that of Section II.

Proof of Theorem II.1.

Let $\Theta,\bar{\Theta}\in\mathcal{D}$. For $0\leq u\leq v\leq m$, we define $\mathcal{N}^{u:v}:\mathbb{R}^{\ell_{u}}\to\mathbb{R}^{\ell_{v}}$ so that $\mathcal{N}^{u:u}(x)=x$ and $\mathcal{N}^{u:v}(x):=\sigma^{(v)}\left(A^{(v-1)}\mathcal{N}^{u:v-1}(x)+b^{(v-1)}\right)$ for $u<v$. Observe that $\mathcal{N}^{u:v}(x)=\mathcal{N}^{c:v}\circ\mathcal{N}^{u:c}(x)$ when $0\leq u\leq c\leq v\leq m$ and that $\mathcal{N}^{\Theta}_{m}(x)=A^{(m)}\mathcal{N}^{0:m}(x)+b^{(m)}$.

For convenience, we replace $\sigma^{(u)}$ by $\tilde{\sigma}$ when no confusion is possible. For all $||v||\leq D$, we set $F(v):=\sum_{i,j}^{n}\mu\left(\left[g(v)g(v)^{T}-\Sigma\right]_{i,j}\right)$. Replicating (32) in Step 5 of Lemma IV.4, we deduce the bound $L_{F}:=2\sqrt{n}D$ for the Lipschitz constant of $F$. Similarly following the proof of Lemma IV.5 with $C=0$, we can deduce one for $\frac{\partial F}{\partial v_{\nu}}$ as $L_{F^{\prime}}:=2n(1+\mu^{\prime\prime}_{\max}D^{2})$ by (25) and Lemma IV.3.

We fix $0\leq u\leq m$, $1\leq i\leq\ell_{u+1}$, and $1\leq j\leq\ell_{u}$. As in Lemma IV.1, we can check that

[TABLE]

Note that $\tilde{\sigma}^{\prime}(v)$ is an $\ell_{u+1}\times\ell_{u+1}$ diagonal matrix. Since $||\mathcal{N}^{0:u}(h(\Sigma))||^{2}\leq\ell_{u}$ for $u>0$, we get

[TABLE]

For $u=0$, we need to replace $\ell_{u}$ by $||\Sigma||^{2}_{S}$. Denoting $\mathcal{N}^{\Theta}_{u}:=A^{(u)}\mathcal{N}^{0:u}(h(\Sigma))+b^{(u)}$ and expanding (37), we have for $u<m$:

[TABLE]

Using exactly the same reasoning as in Step 4 of Lemma IV.4, a Lipschitz constant $\mathcal{L}_{u}$ bound for $\mathcal{N}^{\Theta}_{u}$ satisfies $\mathcal{L}^{2}_{u}:=D^{2}_{+}\mathcal{L}^{2}_{u-1}+\ell_{u}+1$ with $\mathcal{L}^{2}_{0}=1+||\Sigma||_{S}^{2}=L_{Z}^{2}$ and $D_{+}:=\max\{D\sigma^{\prime}_{\max},1.001\}>1$, so that

[TABLE]

To estimate the Lipschitz constant of (39), Lemma IV.3 with $\phi_{1}\equiv\frac{\partial F}{\partial v_{\nu}}$, $\phi_{2}\equiv\mathcal{N}^{\Theta}_{m}$, $\phi_{3}\equiv A^{(m)}_{\nu,i}$ gets us a Lipschitz constant bound of $|A^{(m)}_{\nu,i}|L_{F^{\prime}}\mathcal{L}_{m}$ for $A^{(m)}_{\nu,i}\frac{\partial F(\mathcal{N}^{\Theta}_{m})}{\partial v_{\nu}}$ and a bound of $|A^{(m)}_{\nu,i}|L_{F}$. With $\phi_{1}\equiv\tilde{\sigma}^{\prime}$, $\phi_{2}\equiv\mathcal{N}^{\Theta}_{m-1}$, and $\phi_{3}\equiv A^{(m)T}_{:,i}F^{\prime}(\mathcal{N}^{\Theta}_{m})$, we get a Lipschitz constant estimate of $\left(\sigma^{\prime\prime}_{\max}\mathcal{L}_{m-1}L_{F}+\mathcal{L}_{m}L_{F^{\prime}}\right)||A^{(m)}_{:,i}||$ for $[\sigma^{\prime}(\mathcal{N}^{\Theta}_{m-1})A^{(m)T}F^{\prime}(\mathcal{N}^{\Theta}_{m})]_{i}$. Pursuing as follows, we end with a Lipschitz constant bound for $\frac{\partial\varphi(\Theta)}{\partial b_{i}^{(u)}}$ and $u<m$ of

[TABLE]

For $u=m$, we can simply take $L_{i}^{(m)}:=L_{F^{\prime}}$. Note that

[TABLE]

Using (38) and summing up over $i$ yields

[TABLE]

Finally, by summing $u$ from [math] to $m$, we obtain that

[TABLE]

where $\mathcal{C}$ is a constant that only depends polynomially on $\mu_{\max}^{\prime\prime}$, $\sigma_{\max}^{\prime}$, and $\sigma_{\max}^{\prime\prime}$ (with powers in $O(m)$).

The second part of Theorem II.1 is well-known in optimization theory; see, e.g., [18, Section 1.2.3]. ∎

Acknowledgment

C. Herrera gratefully acknowledges the support from ETH-foundation and from the Swiss National Foundation for the project Mathematical Finance in the light of machine learning, SNF Project 172815.

A. Neufeld gratefully acknowledges the financial support by his Nanyang Assistant Professorship (NAP) Grant Machine Learning based Algorithms in Finance and Insurance.

The authors would like to thank Florian Krach, Hartmut Maennel, Maximilian Nitzschner, and Martin Stefanik for their careful reading and suggestions. Special thanks go to Josef Teichmann for his numerous helpful suggestions, ideas, and discussions.