4
$\begingroup$

I need to compute the pseudo-inverse of a very large rectangular dense matrix without any special structure or properties. I run out of memory/computing power and have no access to a large parallel computing resource.

However, the good news is that I need only one column of the result at a time (for the subsequent calculations).

Is there any iterative algorithm that can compute the 'kth' column (or at-least progressively build up the first 'k' columns of the pseudo-inverse ?). I'd appreciate any inputs/thoughts on this.

PS: I am using MATLAB for now, but the programming environment does not really matter.

Improve this question
$\endgroup$
2
  • $\begingroup$ Is the matrix short and fat or is it tall and thin (that is, does it have many more columns than rows or vice versa)? Is it of full rank in the smaller dimension? $\endgroup$
    – Brian Borchers
    Commented Oct 14, 2015 at 12:16
  • $\begingroup$ You haven't explained why you want columns of the pseudoinverse as opposed to just solving the inverse problem for some collection of right hand sides. $\endgroup$
    – Brian Borchers
    Commented Oct 15, 2015 at 15:55

1 Answer 1

Reset to default
3
$\begingroup$

The pseudo inverse $A^+$ fulfills: $A^+b$ is the minimimum norm solution of the least squares problem $\min_x \|Ax-b\|_2^2$. Hence, to calculate $A^+e_k$ by solving the respective optimization problem (which has a different structure, depending on the format of $A$, as Brian Borchers commented.

You could also use an iterative method to solve $A^TA x = A^Te_k$ in the overdetermined case or solve $AA^Ty = e_k$ and $x = A^T y$ in the underdetermined case, respectively. The method of choice depends on further properties of the matrix…

Improve this answer
$\endgroup$
6
  • 1
    $\begingroup$ If $A$ isn't of full row rank and also not of full column rank, then the normal equations will be singular- you would need to regularize the problem to get a minimum norm least squares solution. $\endgroup$
    – Brian Borchers
    Commented Oct 14, 2015 at 14:46
  • $\begingroup$ LSQR could also be used for the problem - in exact arithmetic it is the same as CG on the normal equations, but for ill conditioned problems in inexact arithmetic it has some prefferable properties. See, web.stanford.edu/class/cme324/paige-saunders2.pdf $\endgroup$
    – Nick Alger
    Commented Oct 14, 2015 at 17:38
  • $\begingroup$ @BrianBorchers Good point! But wouldn't CG, initialized with zero, also converge to the right solution of the normal equations $A^TA x = A^T e_k$? (A solution exists, even if $A^TA$ is not of full rank, since $A^Te_k$ is in the range of $A^TA$.) $\endgroup$
    – Dirk
    Commented Oct 15, 2015 at 7:35
  • $\begingroup$ CG applied to the normal equations has the property that the norm of the solution increases monotonically. This has a useful regularizing effect, but only if you truncate the CG iterations early. On a singular problem like this, if you keep iterating, the norm of the solution will blow up. The LSQR code has an option for explicitly regularizing the solution, and that's what I would probably use. $\endgroup$
    – Brian Borchers
    Commented Oct 15, 2015 at 15:38
  • $\begingroup$ @BrianBorchers Really? If the system is well posed in the sense that the ratio of the largest and the smallest nonzero singular value is not too large, CG does not have problems if the system has a solution. In this problem, the system has a solution, by construction. CG should converge and the number of steps (in exact arithmetic) should be equal to the dimension of the range. In practice you'll not see a difference for well conditioned systems. I am not sure from where I learned this (probably Engl/Hanke/Neubauer?). $\endgroup$
    – Dirk
    Commented Oct 17, 2015 at 18:44

Your Answer

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

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