6
$\begingroup$

I have an $n \times n$ unsymmetric matrix $A$ that results from the discretization of an ill-posed Poisson problem, and thus is rank-deficient with null space of dimension one. I want to compute just the smallest singular triplet

\begin{equation*} (u_n,\sigma_n,v_n) \end{equation*}

where

\begin{equation*} u_n^T A v_n = \sigma_n = 0. \end{equation*}

I have found, through some experimentation, that applying inverse iteration to the shifted (and thus full-rank) matrix $A^T - \sigma I$ yields the left singular vector very quickly, often in one or two iterations. In Matlab code: 

x(:,1) = rand(n,1); 
for k = 1:MAXIT
    x(:,k+1) = ( A' - sigma*eye(n) )\x(:,k); 
    x(:,k+1) = x(:,k+1)/norm(x(:,k+1)); 
end

Usually, this iteration converges to small enough tolerances ($|| x_k^T A|| = 10^{-9}$ or so) in one or two iterations.

I think this converges so rapidly because the largest and second largest singular values of $(A^T - \sigma I)^{-1}$ are spaced by $\sigma^{-1}$ which can be made large, but I have not been able to prove that this inverse iteration should converge. I'm sure there is some clever way of using the SVD to show this convergence as in the standard proof of the equivalent inverse iteration for the eigenvectors, but I can't come up with it.

Any help, ideas, or thoughts on why this works, and how I could prove this works would be helpful.

Improve this question
$\endgroup$
9
  • 1
    $\begingroup$ What do you tend to use for $\sigma$? Inverse iteration should converge to the eigenvector associated with the eigenvalue closest to $\sigma$. Apart from that, if you know $\sigma_n = 0$, then the left singular vector should be the eigenvector associated with the zero eigenvalue of $A^*$. Inverse iteration converges with rate $|\sigma - \lambda_{\text{closest to }\sigma}|/|\sigma-\lambda_{\text{next closest to }\sigma}|$. Can you say anything about the eigenvalues of $A^*$ (or equivalently, $A$)? $\endgroup$
    – Jesse Chan
    Jan 6, 2014 at 15:25
  • $\begingroup$ What method do you use to compute the inverse matrix-vector product in the inverse power iteration? If it is anything pivoted, then applying the same pivoted factorization method to the unmodified matrix $A$ should give you the kernel vectors using the last row and column of the factorization. $\endgroup$
    – Lutz Lehmann
    Jan 6, 2014 at 17:27
  • 2
    $\begingroup$ $v_n$ is a kernel vector of $A$ of norm $1$, the same for $u_n$ and $A^\top$. If you compute $A=PLU$, then $Ax=0$ is equivalent to $Ux=0$ and by assumption and construction, the last row of $U$ is zero, which allows to set $x_{n}=1$ and solve the other equations in that triangular system. Then $v_n=x/\|x\|$. $\endgroup$
    – Lutz Lehmann
    Jan 6, 2014 at 19:49
  • 1
    $\begingroup$ Pivoting in the QR decomposition requires column switches, so $AP=QR$ is the result, with $R$ again upper triangular and last diagonal entry equal to zero. Thus again, solve $Rx=0$ with $x_n=1$ and set $v_n=Px/\|x\|$. The left singular eigenvector is even easier to compute, it is $u_n=y/\|y\|$ with $y=Q^\top e_n$. $\endgroup$
    – Lutz Lehmann
    Jan 6, 2014 at 19:56
  • 1
    $\begingroup$ Sumedh, I meant that if you have a zero singular value, this implies you have a zero eigenvalue as well - if $Ax = 0$, $A^*Ax = 0$ as well. Thus, your inverse iteration on $A^*$ is simply converging to the zero eigenvector, which is why it works. LutzL notes this as well in his PLU comments. $\endgroup$
    – Jesse Chan
    Jan 6, 2014 at 20:28

1 Answer 1

Reset to default
3
$\begingroup$

Just collecting the answers contained in the comments:

The left and right singular vectors of $A$ are exactly the left and right kernel vectors of $A$. Since the assumption is that the kernels have dimension $1$ and $A$ is well conditioned outside the kernel, one can compute the kernel vectors using the usual, pivoted, matrix factorizations.

For the LU factorization, $A=PLU$, $P$ a permutation matrix, $L$ lower triangular with $1$-diagonal, $Ax=0$ is equivalent to $Ux=0$. The last row of $U$ including the diagonal element is zero, thus choosing $x_n=1$ and substituting backwards results in a kernel vector $x$ and $v_n=x/\|x\|$.

$y^TA=0$ is equivalent to $(y^TPL)R=0$, which has the nontrivial solution $(y^TPL)=e_n^T$ or $L^T(P^Ty)=e_n$. $L^T$ is a regular triangular matrix, thus this is easily solvable and $y=P(L^{-T}e_n)$.

In the QR decomposition, row pivoting is without effect, so column pivoting needs to be applied, the factorization is then $AP=QR$. $Ax=0$ is equivalent to $R(P^Tx)=0$ which can be solved as before for $P^Tx$, and $y^TA=0$ is equivalent to $(y^TQ)R=0$, again with the obvious solution $y^TQ=e_n^T$ or $y=Qe_n$

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks for summarizing this for future readers. As a matter of courtesy, is it the role of the question asker to do this? I'm just getting my feet wet in scicomp, and don't know the rules of the road yet :) $\endgroup$
    – Sumedh Joshi
    Jan 7, 2014 at 16:11
  • $\begingroup$ No, it is the task of the responders to decide when a comment turns from requesting corrections or more information on the problem formulation to providing partial or full answers. That is, I should have put my second comment into an answer from the start, but it was not clear to me if there were other reasons to avoid direct factorization. $\endgroup$
    – Lutz Lehmann
    Jan 7, 2014 at 16:43
  • $\begingroup$ I think you need complete pivoting instead of partial pivoting here. Otherwise the $LU$ factorization fails for $A = \begin{pmatrix} 0 & 2 \\ 0 & 1 \end{pmatrix}$. $\endgroup$
    – wim
    Apr 20, 2019 at 8:19
  • $\begingroup$ .... Unless the $LU$ factorization code simply neglects the zero column, which seems to be the case in Octave. $\endgroup$
    – wim
    Apr 22, 2019 at 8:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.