# Probability distribution associated with nuclear norm?

The $$\ell_1$$ norm of model parameters is often added to loss functions because it induces sparsity in the solution of the overall cost function:

$$c(\theta) = -\log L(x|\theta) + \lambda ||\theta||_1$$

The fact that the added loss term can be written as the log of the density of a Laplace or Double Exponential distribution reveals a Bayesian perspective on this procedure. The Laplace density $$\delta_L$$ is given by:

$$\delta_L \propto e^{-\lambda ||\theta||_1}$$

Similarly, the nuclear norm $$||A||_*$$ is increasingly being used in cost functions due to the low rank property it induces in solutions. I would like to determine what the analagous Bayesian procedure would be, in terms of what prior it would place on the unknown $$A$$.

The matrix distribution would thus be of the form $$\delta_N \propto e^{-\lambda||A||_*}$$. Is this a well studied distribution?

• Did you find anything else on this? Since the nuc norm is the sum of single values of A, you get a product of exponential distributions on those. Implies the normalization constant has a term $\lambda^{-\min(n,m)}$ when A is n x m? Commented Nov 3, 2021 at 16:46
• @daknowles Unfortunately I still don't have much information about this distribution. Commented Nov 3, 2021 at 18:49

We can derive the normalizing constant $$C(\lambda):=\int e^{-\lambda ||M||_*}dM$$ using the smooth coarea formula.

Suppose we are given smooth functions $$F: \mathbb{R}^n\to\mathbb{R}$$ and $$f: \mathbb{R}\to\mathbb{R}$$. The formula states that $$\int_{\mathbb{R}^n} f(F(x))dx=\int_{t\in \mathbb{R}}\left(\int_{x\in F^{-1}t} f(x) |\nabla F(x)|^{-1}d H_t(x) \right)dt$$

where $$dH_t$$ denotes the surface measure on the preimage $$F^{-1}t$$. In our case, we'll have $$f(x)=e^{-\lambda x}$$ and $$F(M)=||M||_*$$. Also, denote the dimensions of $$M$$ by $$n,m$$

So first of all we'll need to compute $$|\nabla F(M)|$$. By the resuls from these notes, we can write down the derivatives of the singular values of a matrix. Suppose we have an SVD decomposition $$M=USV^t$$, with $$U^tU=I_{min(n,m)}=V^tV$$ and $$S$$ the diagonal matrix of singular values (we can assume $$M$$ is full rank, because the corresponding set of matrices has full measure).Then the differential with respect to the ith singular value is given by $$dS_i=(U^tdMV)_{ii}$$. Summing over these differentials, and recalling that the nuclear norm is the sum of singular values, we get $$\begin{eqnarray*} d||M||_{*}&=&Tr(U^tdMV)\\ &=&\sum_i \sum_{a,b}U[a,i]V[b,i]dM_{ab}\\ & = & \sum_{a,b} dM_{ab} \sum_i U[a,i]V[b,i]\\ & = & \sum_{ab}dM_{ab} (UV^t)_{ab} \end{eqnarray*}$$

In other words, if we consider the gradient $$\nabla ||M||_*$$ with respect to the entries of $$M$$, then the $$ab$$ coefficient is given by $$(UV^t)_{ab}$$. Now computing the norm of the gradient is straight forward:

$$\begin{eqnarray*} |\nabla ||M||_*|^2& = & \sum_{ab} (UV^t)_{ab}^2\\ & = & Tr((UV^t)(UV^t)^t)\\ & = & Tr(UV^tVU^t)\\ & = & Tr(U^tUV^tV)\\ & = & Tr(I_{min(n,m)}I_{min(n,m)})\\ & = & min(n,m) \end{eqnarray*}$$ In the second equality we used the fact that $$\sum_{ab} X_{ab}^2=Tr(XX^t)$$ for any matrix $$X$$.

Now, we can compute $$\int e^{-\lambda F(M)}dM$$ using the coarea formula (and recall that $$F(M):=||M||_*$$). We get $$\int_0^{\infty} \left(\int_{M\in B^*(t)} |\nabla F(M)|^{-1}e^{-\lambda F(M)} dM_{|B^*(t)}\right)dt$$

where $$B^*(t)$$ is the set of all matrices of nuclear norm $$t$$, and $$dM_{|B^*(t)}$$ denotes the restriction of the Euclidean measure to this set. Conveniently, the integrand is constant over each preimage, $$min(n,m)^{-1/2}\int_0^{\infty} e^{-\lambda t} \int_{B^*(t)} dM_{|B^*(t)}dt=min(n,m)^{-1/2}\int_0^{\infty} e^{-\lambda t} Vol(B^*(t))$$

The ball $$B^*(t)$$ is a dilation of the unit ball $$B^*(1)$$,so the volume is $$Vol(B^*(1)) t^{nm-1}$$.

Finally, we get $$min(n,m)^{-1/2}Vol(B^*(1))\int_0^{\infty}e^{-\lambda t}t^{nm-1}dt=min(n,m)^{-1/2}Vol(B^*(1))(nm-1)!\lambda^{-nm}$$

for the normalizing constant $$C(\lambda)$$.

Given this, we can also derive the distribution of the random variable $$||M||_*$$. To start with, we'll take the identity $$\int e^{-\lambda ||M||_*}=C(\lambda)$$ and differentiate under the integral to see:

$$(-1 )^k\int ||M||^k_* e^{-\lambda ||M||_*}=d^k_{\lambda}C(\lambda)=(-1)^k{\frac {(nm)^{(k)}}{\lambda^k}}C(\lambda)$$

where $$(nm)^{(k)}={\frac {(nm+k-1)!}{(nm-1)!}}$$ denotes the rising factorial. Thus $$E ||M||^k_*=C(\lambda)^{-1}\int ||M||^k_* e^{-\lambda ||M||_*}dM= \lambda^{-k}{\frac {(nm+k-1)!}{(nm-1)!}}=\lambda^{-k}{\frac {\Gamma(nm+k)}{\Gamma(nm)}}$$

Remarkably, these are the same moments as a $$\Gamma$$ distribution with rate parameter $$\lambda$$ and shape parameter $$mn$$. By uniqueness of moment generating functions, we conclude that $$||M||_*\sim \Gamma(mn,\lambda)$$.

You can also use a similar argument to derive the distribution over singular vectors/values. To phrase the question, suppose that I have sampled orthogonal matrices $$U,V$$ independently and uniformly at random, as well as a positive vector $$\sigma$$ from some fixed distribution $$p(x)$$ on $$\mathbb{R}_+^n$$. The question is: What should $$p$$ be, in order to ensure that the resulting matrix $$M:=UDiag(\sigma)V^t$$ is distributed according to the original distribution $$e^{-\lambda ||M||_*}$$?

Assume WLOG that $$m>n$$. By the change of variables formula, the resulting density of $$M$$ is given by $$|{\frac {\partial (U,V,\sigma)}{\partial M}}|p(\sigma)$$,were $$U,V,\sigma$$ are the SVD of $$M$$.

Lemma $$\left|{\frac {\partial (U,V,\sigma)}{\partial M}}\right|=\prod_{i<=n} \sigma_i^{-(m-n)}\prod_{i

proof

See Aspects of Multivariate Statistical Theory by R.J. Muirhead. I also worked it out in a previous edit to this question, so it is visible in the edit history. It is not particularly difficult, provided that one uses the formulas for the derivative of the SVD in the linked notes. There are two tricks that make it considerably easier: (1) by symmetry it suffices to assume that $$U=I$$ and $$V=\left (\begin{array}{c}I\\ 0\end{array}\right)$$, (2) rather than the raw jacobian, it turns out to be easier to work with the Gram matrix $$G_{ab}=\sum_{ij} {\frac {\partial x_a}{\partial M_{ij}}}{\frac {\partial x_b}{\partial M_{ij}}}$$ where $$x_a$$ ranges over all independent entries of $$U,V,\sigma$$. This matrix ends up having a particularly simple structure.

So in summary:

The original distribution $$e^{-\lambda ||M||_*}$$ is obtained as the distribution of a matrix of the form $$UDiag(\sigma)V^T$$, where $$U$$ and $$V$$ are sampled uniformly from the space of all orthogonal matrices, and $$\sigma$$ is sampled from $$\mathbb{R}_+^n$$ according to the density $$\propto \prod_{i\leq n} \sigma_i^{m-n}\prod_{i.

singular value distribution

To explore the singular values, I generated samples from the singular value distribution using MCMC. in this case, $$n=m=20$$ and $$\lambda=1$$. Below is the histogram of the resulting sums $$\sum_i \sigma_i$$. These sums appear to follow a $$\Gamma(20*20,1)$$ distribution (shown in orange), as predicted by the above discussion.

A bit more interesting is the distribution over the individual singular values. If $$\sigma_k$$ denotes the $$k$$th smallest singular value, then it appears that each $$\sigma_k$$ is in fact well approximated by a gamma distribution. Below are plotted histograms for a range of $$\sigma_k$$, with the approximating gamma pdfs overlaid.

This is quite interesting, and I am not sure why it would be the case.

analysis of largest singular value

It turns out that there is a relatively simple analytic expression for the largest singular value, which can be found by applying the techniques from this paper. The key property is that the singular value density looks like $$\propto |det \{\phi_j(\sigma_i)\}_{ij}|$$ where $$\phi_i$$ are scalar-valued functions. Due to the Vandermonde determinant formula, it is easy to see that we can take $$\phi_j(\sigma)=\sigma^{m-n}e^{-\lambda\sigma}\sigma^{2(j-1)}$$.

To evaluate the cdf of the largest singular value, we need to evaluate the integral $$\int_{0<\sigma_1<\dots<\sigma_n. In general, it turns out that an integral of this form can be expressed as $$\sqrt{|det(A)|}$$ where $$A_{ij}=\int_{[0,t]^2}sgn(x-y)\phi_i(x)\phi_j(t)dxdy$$ (for even $$n$$; if n is odd the definition of $$A$$ is slightly more complicated). Even better, it turns out that, for our $$\phi$$ functions, the entries of $$A$$ satisfy a recurrence relation that obviates the need to actually evaluate the integral. See the linked paper for more details.

• Thanks for all the hard work. Unfortunately, my smooth manifolds theory is quite limited, so sorry if this is a dumb question, but our function $F$ is not smooth, right? But you're saying something along the lines of "since it is nonsmooth over a set of eucliden measure zero, we can treat it as effectively smooth"? In any event, the implied distribution over the singular values is very helpful, thanks. So interesting that they are dependent. Commented Jan 7, 2022 at 18:37
• @John Madden, Yes you're right, if you restrict to the set of matrices with no repeated singular values, then $F$ will be smooth (since for each $k$, the function $M\mapsto \sigma_k(M)$ is well-defined and smooth on that set). Like you said, that's good enough because you are only throwing away a measure zero set. I'm actually not sure what happens on that set, as far as the smoothness of $F$. Commented Jan 7, 2022 at 19:26