8
$\begingroup$

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.

sparse Jacobian

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – R zu Oct 25 '18 at 15:10
  • 3
    $\begingroup$ 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. $\endgroup$ – GoHokies Oct 25 '18 at 16:25
  • 1
    $\begingroup$ 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. $\endgroup$ – GoHokies Oct 25 '18 at 16:30
  • 2
    $\begingroup$ @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.) $\endgroup$ – Christian Clason Oct 25 '18 at 16:35
  • 2
    $\begingroup$ @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! $\endgroup$ – Richard Zhang Oct 26 '18 at 16:52

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Browse other questions tagged or ask your own question.