# 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