# How to find a good preconditioner to the system $(A^T A + \lambda I) x = A^T b$?

The system in the title has a damper factor $$\lambda > 0$$ and the matrix $$A$$ is sparse and rectangular, with a structure I can exploit to solve matrix vector products very fast. My current solver, LSMR, is trying to solve the normal equations $$(A^TA + \lambda I) x = A^T b$$ associated to the original problem $$\min \|Ax - b\|$$.

Although each iteration is computed very fast, the algorithm uses the maximum number of iterations. I know this can be fixed with a good preconditioner. This is where lies my problem.

$$A^TA + \lambda I$$ is a SPD matrix, which is a good property to have. On the other side, this matrix is no more sparse. I don't know how to choose and use a preconditioner for this dense matrix. I suppose this is already worked by someone.

I want to know how to proceed in this case and, if possible, how to use the sparsity of $$A$$ to obtain a good preconditioner. What are the common approaches?

EDIT: In order to be more complete, I'll briefly describe how the matrix $$A$$ is obtained. My problem at hand consists in minimizing the error associated to a low tensor rank-$$r$$ approximation. You can consider a tensor $$T$$ as being a multidimensional array. In this case, a 3-D multidimensional array with coordinates $$T_{ijk}$$, for $$1 \leq i \leq m, 1 \leq j \leq n, 1 \leq k \leq p$$. I am considering an approximation $$\tilde{T}_{ijk} = \sum_{\ell=1}^r X_{i \ell} \cdot Y_{j \ell} \cdot Z_{k \ell}$$. The error in this approximation is given by $$\frac{1}{2} \sum_{i,j,k} \left( T_{ijk} - \tilde{T}_{ijk} \right)^2 = \frac{1}{2} \sum_{i,j,k} res_{i,j,k} (X,Y,Z)^2,$$ where $$X, Y, Z$$ lists all components $$X_{i \ell}, Y_{j \ell}, Z_{k \ell}$$ and $$res_{ijk}$$ is the residual of the component with index $$i,j,k$$.

To find the components of $$\tilde{T}$$ which minimize the error above, it is of interest to find the Jacobian matrix of $$res = (res_{111}, res_{112}, \ldots, res_{mnp})$$. We have the formulas below for the partial derivatives:

$$\frac{\partial res_{ijk}}{\partial X_{I \ell}} = \left\{ \begin{array}{c} - Y_{j \ell} Z_{k \ell},\quad \text{if } i = I,\\ 0, \quad \text{otherwise} \end{array}\right.$$

$$\frac{\partial res_{ijk}}{\partial Y_{J \ell}} = \left\{ \begin{array}{c} - X_{j \ell} Z_{k \ell},\quad \text{if } j = J,\\ 0, \quad \text{otherwise} \end{array}\right.$$

$$\frac{\partial res_{ijk}}{\partial Z_{K \ell}} = \left\{ \begin{array}{c} - X_{i \ell} Y_{j \ell},\quad \text{if } k = K,\\ 0, \quad \text{otherwise} \end{array}\right.$$

This will give a sparse matrix, which becomes more sparse as we increase the dimensions. The structure follows a nested for loop pattern, from left to right. The figure below shows this structure for $$m = 3, n = 5, p = 7, r = 10$$. I hope this can be useful for someone to spot the "right" preconditioner, because at the moment I really don't know how to proceed. Keep im mind that I'm trying to use this structure to find a preconditioner for $$A^TA + \lambda I$$, where $$A$$ is this sparse matrix just described.

• Normal equation squares the condition number. Other methods are better if the matrix is ill-condition. (see: eigen.tuxfamily.org/dox/group__LeastSquares.html) I haven't studied how to combine preconditioning with the other methods though. – R zu Oct 25 '18 at 15:10
• have you tried the "usual suspects", i.e. diagonal preconditioning or incomplete Cholesky? other, less known preconditioning strategies for sparse LS are evaluated in this recent ACM TOMS paper by Gould. – GoHokies Oct 25 '18 at 16:25
• also, it would be good to include more information about the origin and (block-)structure of $A$. good preconditioners are often built from such domain knowledge. – GoHokies Oct 25 '18 at 16:30
• @Rzu Note that Integral mentioned they were using LSMR, not CG for the normal equation, which avoids the numerical issues with the squared condition number. (@Integral: LSMR is equivalent to MINRES applied to the normal equations, not CG -- that's a different method.) – Christian Clason Oct 25 '18 at 16:35
• @Integral As someone who works in this field, I would say this is basically an open research question. With a good preconditioner for the trust-region subproblem of an outer tensor decomposition, you will give all the SGD people a run for their money. Good luck! – Richard Zhang Oct 26 '18 at 16:52

Have you tried the preconditioner diag(A'*A)+lmb*speye(size(A))? For lmb large enough, I would imagine this preconditioner would accelerate the convergence significantly.
diag(A'*A) can be computed fairly efficiently using $$B(i,i) = \sum_j A^T(i,j)A(j,i) = \sum_j A(j,i)A(j,i)$$.
Another idea is to use fixed point iterations as preconditioners. Since your matrix is SPD, 4 steps of Gauss-Seidel should be a decent preconditioner. However, avoiding the explicit construction of the lower triangular part of $$A^TA$$ is challenging. It is doable though and I know that few of my old colleagues had implemented those for personal use.
• I'm not familiar with the specifics of LSMR, but depending on how and where the preconditioner is applied, it may be necessary to preserve symmetry of the matrix. You can adapt this preconditioner to that caes by muliplying from the left and the right by $B^{1/2}+\sqrt{\lambda}I$, where $B$ is as constructed above. This will most likely have a similar effect as a preconditioner while preserving symmetry – whpowell96 Feb 16 at 1:39