2
$\begingroup$

Let us assume that we have a function, $f(A)=\text{vec}(A^{-1})^\intercal B$, dependent on $A^{-1}$. However, due to some machine-precision limitations, the programming language I'm using cannot invert $A$ even if it has $Det(A)>0$. So, I'm thinking of using a pseudo-inverse, instead.

However, I don't know if there's a 'distance' relationship between the true inverse and the pseudo-inverse. I just know that when there's an inverse, both are equal. But it would be nice, that when the original matrix A and the 'perturbed/truncated' version of A are close, that the pseudo-inverse be close to the inverse of the original...

I'm also looking at pseudo-inverses for which there's a relatively easy and efficient way to compute.

Any suggestions?

Edit: $\dim(A)=DT\times DT$ and $\dim(B)=D^2T^2\times D^2$

Following Federico Poloni's comment, I think I could actually recompute the function above as $$\left[\cdots Tr\left(\text{LinearSolve}\left(A^\intercal, C_i \right) \right) \cdots \right]_{1\times i=1,...,D^2}$$ where $C_i=\text{Reshape}_{DT\times DT}\left(\text{Column}(B,i)\right)$

$\endgroup$
9
  • $\begingroup$ what's the action of $\rm{vec}(\cdot)$? is it to reshape the inverse into a one-dimensional (column?) vector? $\endgroup$
    – GoHokies
    Commented Apr 27, 2019 at 8:14
  • $\begingroup$ @GoHokies Yes. It's an example of a function, among others, which I use with reshapping operations. ;) $\endgroup$ Commented Apr 27, 2019 at 8:23
  • 1
    $\begingroup$ ok, so if $A$ is $N \times N$, then $B$ is an $N^2 \times N^2$ matrix - or a $N^2 \times 1$ vector? $\endgroup$
    – GoHokies
    Commented Apr 27, 2019 at 8:27
  • 1
    $\begingroup$ Is $\det(A)>0$ coming from an exact formula, or is it just a numerical observation? Determinants give no easily readable information on the distance of a matrix from singularity, so in numerical practice they are seldom used. $\endgroup$ Commented Apr 27, 2019 at 8:37
  • 1
    $\begingroup$ Also: you write "no linear solve", but what you have here can be rewritten as linear solves: vectorize (each column of) $B$ to obtain a $N\times N$ matrix $C$, then compute $A^{-1}C$, and the quantity you are looking for is its trace. $\endgroup$ Commented Apr 27, 2019 at 18:37

1 Answer 1

6
$\begingroup$

Pseudoinverses typically will be computed via some truncation procedure to determine the rank, so they are not close to the original inverse. Example: $$A = Q \begin{bmatrix} 1\\ & 10^{-12} \\ & & 10^{-18} \end{bmatrix}Q^{*} $$ has exact (pseudo)inverse $$A^+ = A^{-1} = Q \begin{bmatrix} 1\\ & 10^{12} \\ & & 10^{18} \end{bmatrix}Q^{*}, $$ but any numerical procedure will have to take a decision on the rank and will likely truncate one or two of those diagonal entries to zero, returning for instance $$B = Q \begin{bmatrix} 1\\ & 10^{12} \\ & & 0 \end{bmatrix}Q^{*}, $$ or $$B = Q \begin{bmatrix} 1\\ & 0 \\ & & 0 \end{bmatrix}Q^{*}, $$ (depending on truncation thresholds), which are not close at all to $A^+$. That is an inherent limitation and I don't think you can overcome it. It's by design: typically, when one writes pinv(A) one wants $B$, not $A^{-1}$. If you really want $A^{-1}$, then just use inv(A) and forget about the warnings.

But it's probably time to go back to the blackboard, think about what you really need in these edge cases, and see if there is a more stable alternative expression for it. If the exact solution of your problem is that ill-conditioned, then probably it's not even what you want at the end of the day.

$\endgroup$

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.