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 agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.