An Efficient Approach for Computing Optimal Low-Rank Regularized Inverse Matrices

By interpreting problem (3) as an empirical Bayes risk minimization problem [18, 13, 19], we can provide a statistical framework for analyzing the problem and motivate numerical algorithms for computing an optimal low-rank regularized inverse matrix. In the following we assume that ${\boldsymbol{\xi}}$ and ${\boldsymbol{\delta}}$ are random variables with given probability distributions $\mathcal{P}_{{\boldsymbol{\xi}}}$ and $\mathcal{P}_{{\boldsymbol{\delta}}}$. Using the fact that ${\bf Z}{\bf b}-{\boldsymbol{\xi}}=({\bf Z}{\bf A}-{\bf I}){\boldsymbol{\xi}}+{\bf Z}{\boldsymbol{\delta}}$, the Bayes risk minimization problem can be written as,

where $\mathbb{E}_{{\boldsymbol{\xi}},{\boldsymbol{\delta}}}$ is the expected value with respect to the joint distribution of ${\boldsymbol{\xi}}$ and ${\boldsymbol{\delta}}$ [18, 13]. The particular choices of $\rho,$ $\mathcal{P}_{{\boldsymbol{\xi}}},$ and $\mathcal{P}_{{\boldsymbol{\delta}}}$ in (5) determine an optimal reconstruction matrix $\check{\bf Z}$. Optimization problem (3) can be interpreted as an empirical Bayes risk minimization problem corresponding to Bayes risk minimization problem (5), where ${\boldsymbol{\xi}}^{(k)}$ and ${\bf b}^{(k)}={\bf A}{\boldsymbol{\xi}}^{(k)}+{\boldsymbol{\delta}}^{(k)}$ denote realizations of the random variables ${\boldsymbol{\xi}}$ and ${\boldsymbol{\delta}}$.

For the special case where $\rho({\bf x})=\left\|{\bf x}\right\|_{2}^{2},$ the Bayes risk can be written as

In this section, we provide closed form solutions for both the unconstrained and rank-constrained Bayes risk minimization problems for the 2-norm case.

We assume that the random variables ${\boldsymbol{\xi}}$ and ${\boldsymbol{\delta}}$ are statistically independent. Let ${\boldsymbol{\mu}}_{{\boldsymbol{\xi}}}$ be the mean and ${\bf C}_{{\boldsymbol{\xi}}}$ be the covariance matrix for the distribution $\mathcal{P}_{{\boldsymbol{\xi}}}$. Similarly, ${\boldsymbol{\mu}}_{{\boldsymbol{\delta}}}$ and ${\bf C}_{{\boldsymbol{\delta}}}$ are the mean and covariance matrix for $\mathcal{P}_{{\boldsymbol{\delta}}}$. We assume that both matrices are positive definite and that all entries are finite. One fact that will be useful to us is that the expectation of the quadratic form ${\boldsymbol{\phi}}^{\top}{\bf K}{\boldsymbol{\phi}}$ involving a symmetric matrix ${\bf K}$ and a random variable ${\boldsymbol{\phi}}$ with covariance matrix ${\bf C}_{{\boldsymbol{\phi}}}$ can be written as

where ${\rm tr\!}\left(\,\cdot\,\right)$ denotes the trace of a matrix [20]. With this framework and the Cholesky decompositions ${\bf C}_{{\boldsymbol{\xi}}}={\bf M}_{{\boldsymbol{\xi}}}{\bf M}_{{\boldsymbol{\xi}}}^{\top}$ and ${\bf C}_{{\boldsymbol{\delta}}}={\bf M}_{{\boldsymbol{\delta}}}{\bf M}_{{\boldsymbol{\delta}}}^{\top}$, we can write the Bayes risk (6) as,

where $\left\|\,\cdot\,\right\|_{\rm F}$ denotes the Frobenius norm of a matrix.

We can simplify this final expression by pre-whitening the noise ${\boldsymbol{\delta}}$, thereby converting to a coordinate system where the noise mean is zero and the noise covariance is a multiple of the identity. Define the random variable

Then, ${\boldsymbol{\mu}}_{{\boldsymbol{\chi}}}={\bf 0}$ and ${\bf C}_{{\boldsymbol{\chi}}}={\bf I}_{m},$ and we have

which can be written as,

where $\widetilde{\bf A}={\bf M}_{{\boldsymbol{\delta}}}^{-1}{\bf A}$ and $\widetilde{\bf b}={\bf M}_{{\boldsymbol{\delta}}}^{-1}({\bf b}-{\boldsymbol{\mu}}_{{\boldsymbol{\delta}}}).$ Thus, without loss of generality (although with the need for numerical care if ${\bf M}_{{\boldsymbol{\delta}}}$ is ill-conditioned), we can always assume that ${\boldsymbol{\mu}}_{\delta}={\bf 0}$ and ${\bf C}_{{\boldsymbol{\delta}}}=\eta^{2}{\bf I}_{m}$. Therefore, our Bayes risk (6) takes the form

Under certain conditions, summarized in the following theorem, there exists a unique global minimizer for (10).

To summarize, we have minimized the Bayes risk (6) under the assumptions that ${\boldsymbol{\xi}}$ and ${\boldsymbol{\delta}}$ are statistically independent random variables with symmetric positive definite finite covariance matrices.

Next we consider the rank-constrained Bayes risk minimization problem. With the additional assumption that ${\boldsymbol{\mu}}_{{\boldsymbol{\xi}}}={\bf 0}$, problem (5) can be written as

where ${\bf M}={\bf M}_{{\boldsymbol{\xi}}}$ and $\alpha=\eta.$ A closed form solution to (15) was recently provided in [21], and the result is summarized in the following theorem.

Next we consider the case where ${\bf C}_{{\boldsymbol{\xi}}}=\beta^{2}{\bf I}_{n}$ for some scalar $\beta$. In this case the Bayes risk simplifies to

and, by rescaling, the low-rank Bayes risk minimization problem can be restated as optimization problem (15) where $\alpha^{2}=\frac{\eta^{2}}{\beta^{2}}$ represents the noise-to-signal ratio and ${\bf M}={\bf I}_{n}$. It was shown in [21] that for this problem, the global minimizer simplifies to

where ${\boldsymbol{\Psi}}_{r}={\rm diag\!}\left(\frac{\sigma_{1}}{\sigma_{1}^{2}+\alpha^{2}},\ldots,\frac{\sigma_{J}}{\sigma_{r}^{2}+\alpha^{2}}\right)$ and ${\bf V}_{r}$ and ${\bf U}_{r}$ were defined in the beginning of this section. Moreover, $\check{\bf Z}$ is the unique global minimizer if and only if $\sigma_{r}>\sigma_{r+1}$. For problems where ${\bf A}$ is known, matrix $\check{\bf Z}$ can be used to compute a filtered solution to inverse problems. That is,

where the filter factors are given by

Here, $\check{\bf Z}$ is the truncated Tikhonov reconstruction matrix [21].

For the case $\rho=\left\|\,\cdot\,\right\|_{2}^{2}$, we can use the relationship between the 2-norm of a vector and the Frobenius norm of a matrix to rewrite optimization problem (3) as

where ${\bf B}=[{\bf b}^{(1)},\ldots,{\bf b}^{(K)}]$ and ${\bf C}=[{\boldsymbol{\xi}}^{(1)},\ldots,{\boldsymbol{\xi}}^{(K)}]$. A closed form solution for (19) exists, and the result is summarized in the following theorem.

The main computational cost of constructing $\widehat{\bf Z}$ using Theorem 3.1 lies in two SVD computations (for matrices ${\bf B}$ and ${\bf P}$). If these computations are possible and minimizing the mean squared error is desired, then we propose to use the above result to compute optimal low-rank regularized inverse matrix $\widehat{\bf Z}$, which in turn can be used to solve inverse problems [24].

For general error measures or for large problems with large training sets, we propose a rank-update approach for obtaining optimal low-rank regularized inverse matrix $\widehat{\bf Z}$. The approach is motivated by the following corollary to Theorem 2.2.

Corollary Corollary states that in the Bayes framework for the 2-norm case, a rank-update approach (in exact arithmetic) can be used to compute an optimal low-rank regularized matrix. That is, assume we are given the rank-$J$ approximation $\check{\bf Z}^{J}$ and we are interested in finding $\check{\bf Z}^{J+\ell}$. To calculate the optimal rank-$(J+\ell)$ approximation $\check{\bf Z}^{J+\ell}$, we just need to solve a rank-$\ell$ optimization problem of the form (20) and then update the solution, $\check{\bf Z}^{J+\ell}=\check{\bf Z}^{J}+\bar{\bf Z}^{\ell}$.

We are interested in the empirical Bayes problem, and we extend the rank-update approach to problem (3), where the basic idea is to iteratively build up matrices $\widehat{\bf X}\in\mathbb{R}^{n\times r}$ and $\widehat{\bf Y}\in\mathbb{R}^{m\times r}$ such that $\widehat{\bf Z}=\widehat{\bf X}\widehat{\bf Y}^{\top}.$ Suppose we are given the optimal rank-$J$ matrix $\widehat{\bf Z}$. That is, $\widehat{\bf Z}$ solves (3) for $r=J$. Then, we compute a rank-$\ell$ update to $\widehat{\bf Z}$ by solving the following rank-$\ell$ optimization problem,

and a rank-$(J+\ell)$ regularized inverse matrix can be computed as $\widehat{\bf Z}=\widehat{\bf Z}+\widehat{\bf X}^{\ell}(\widehat{\bf Y}^{\ell})^{\top}$. This process continues until the desired rank is achieved. A summary of the general algorithm can be found below, see Algorithm 1.

For the 1D examples, we use random time signals by generating piecewise constant functions as seen in Figure 2. Here, the number of jumps in each signal is a randomly chosen integer between 3 and 20, drawn from a uniform distribution. The times for the jumps and the signal strength were both randomly selected in the interval $[0,1]$, also from a uniform distribution. The discretized signal has length $n=150.$ The corresponding blurred signals were computed as in Equation (1), where matrix ${\bf A}$ is $150\times 150$ and represents convolution with a 1D Gaussian kernel with zero mean and variance 2, and ${\boldsymbol{\delta}}$ represents Gaussian white noise. The noise level for the blurred signals was $0.01$, which means that $\left\|{\boldsymbol{\delta}}\right\|_{2}^{2}/\left\|{\bf A}{\boldsymbol{\xi}}\right\|_{2}^{2}\approx 0.01$. Training and validation signals were generated in the same manner. Figure 2 displays four sample signals from the validation set and their corresponding observed signals.

Next, we describe a problem from image deblurring, where the underlying model can be represented as (1) where matrix ${\bf A}$ represents a blurring process, and ${\boldsymbol{\xi}}$ and ${\bf b}$ represent the true and observed images respectively. Entries in matrix ${\bf A}$ depend on the point spread function (PSF). For image deblurring examples, matrix ${\bf A}$ is typically very large and sparse. However, under assumptions of spatial invariance and assumed boundary conditions and/or properties of the PSF, structure in ${\bf A}$ can be exploited and fast transform algorithms can be used to efficiently diagonalize ${\bf A}$ [5]. For example, if the PSF satisfies a double symmetry condition and reflexive boundary conditions on the image are assumed, then ${\bf A}$ can be diagonalized using the 2D discrete cosine transform (DCT). Of course, these approaches rely on knowledge of the exact PSF.

Blind deconvolution methods have been developed for problems where the PSF is not known or partially known [32]; however, most of these approaches require strict assumptions on the PSF such as spatial invariance or parameterization of the blur function, or they require a good initial estimate of the PSF. Our approach uses training data and does not require knowledge of ${\bf A}$ and, hence, knowledge of the PSF. Furthermore, we do not require any assumptions regarding spatial invariance of the blur, boundary conditions, or special properties of the blur (e.g., separability, symmetry, etc.). We refer the reader to [24] for comparisons to blind deconvolution methods.