5

How can I efficiently find a matrix $P_A$ that projects onto the kernel of another matrix $A$? That is, given $Ax=0$, is there an "efficient" (without using inverses) way to find $P_A$ such that $P_A x \in ker(A)$, i.e., $A(P_A y) = 0$ for $Ay \neq 0$?

I know that I can do this with the orthogonal projection $Q_A := I - A^{\dagger}A$ for $A^{\dagger} = A^T(AA^T)^{-1}$ when $A$ has full row rank or $A^{\dagger} = (A^TA)^{-1}A^T$ when $A$ has full column rank, but I'm interested in other projections onto $ker(A)$ where we're not required to compute the pseudoinverse because of two reasons

  1. SVD of $A^TA$ or $AA^T$ might be too expensive or unreliable (because of size or poorly conditioning, respectively)
  2. $A$ might be rank deficient (not have full row rank or full column rank)

In these cases, what's the best way to project onto the kernel of $A$? I imagine that we must rely on knowing something about the structure of $A$.

jjjjjj
  • 2,671
  • 1
    How many rows do $M$ and $X$ have? Any hypothesis on their rank? –  Oct 08 '17 at 17:09
  • Roughly, I've got $M \in \mathbb{R}^{100 \times 10,000}$ with full row rank and $X \in \mathbb{R}^{5,000 \times 10,000}$ with full row rank as well. However, poor conditioning on $X^T X$ (~1e20). – jjjjjj Oct 08 '17 at 17:15
  • 1
    Let us agree on the terms : by $\mathbb{S}^{n\times n}$ do you mean the set of ${n\times n}$ symmetrical matrices ? If so, I am puzzled by your text $(I-M^T M)(X^T X + I) \in \mathbb{S}^{n\times n}$, because the product of two symmetrical matrices, here $(I-M^T M)$ and $(X^T X + I) $ is in general not symmetrical. – Jean Marie Oct 08 '17 at 17:20
  • Sorry, you are right, let me correct, thanks. – jjjjjj Oct 08 '17 at 17:22
  • 1
    I would guess $X$ is poorly conditioned to start with and that forming $X^TX$ makes it much worse. Not trying to slack off but is there no alternate approach than to degrade your problem in such a manner? –  Oct 08 '17 at 17:24
  • @LucJ.Bourhis: that's correct about $X$ in my case. I haven't been able to see one yet...would you be able to discuss more specifically in chat? – jjjjjj Oct 08 '17 at 17:28
  • I'd be happy to! But in my time zone (Paris), I need to call it a day I am afraid! –  Oct 08 '17 at 17:32
  • @LucJ.Bourhis: no problem! Would there be another time to possibly coordinate? – jjjjjj Oct 08 '17 at 17:34
  • Tomorrow about noon Paris time? –  Oct 08 '17 at 18:19
  • @LucJ.Bourhis: that'd be great, thanks! Look for you then – jjjjjj Oct 08 '17 at 18:24
  • Let us continue this discussion in chat. – jjjjjj Oct 09 '17 at 10:07
  • Did you find a satisfactory answer? – a06e Sep 09 '21 at 22:31
  • 1
    $d=A^{\dagger}b$ can be computed by solving $Ad=b$ by computing the minimum-norm least squares solution. In julia programming language, if $h=A^{\dagger}b$, then $b=A \setminus b.$ Hence, in julia language, finding $X$ such that $X=A^{\dagger}A$ is done by calculating $A \setminus A$. Up to my knowledge, there is no method which beats this precision. – R. W. Prado Aug 15 '23 at 17:22

0 Answers0