4
$\begingroup$

I wonder if there is a fast algorithm, say ($\mathcal O(n^3)$) for computing the cofactor matrix (or conjugate matrix) of an $N\times N$ square matrix. And yes, one could first compute its determinant and inverse separately and then multiply them together. But how about this square matrix is non-invertible?

I am curious about the accepted answer on StackOverflow.

What would it mean by "This probably means that also for non-invertible matrixes, there is some clever way to calculate the cofactor (i.e., not use the mathematical formula that you use above, but some other equivalent definition)."?

$\endgroup$
  • 3
    $\begingroup$ welcome to scicomp! I'm curious, too: why do you want to numerically compute the co-factors of a singular matrix (i.e. what's your target application)? $\endgroup$ – GoHokies Jul 9 at 11:58
  • $\begingroup$ I am implementing an algorithm on a large-scale setting that computes conjugate matrix in iterative steps. And I find it is the bottleneck. Maybe there is some mathematical foundation for why it has to be conjugate (instead of others, like pseudoinverse * pseudo-determinant, which all led to divergence on small-scale instance). I am just curious if some one have developed an efficient way to deal with this computational issue. $\endgroup$ – ZUN LI Jul 9 at 14:36
  • 1
    $\begingroup$ If I understand the goal, you want to compute all $n^2$ cofactors with $O (n^3)$ effort. The structure of the underlying matrix will be needed to evaluate the stability of such algorithms. $\endgroup$ – hardmath Jul 10 at 14:16
3
$\begingroup$

So, a cofactor matrix is a transpose of an adjugate matrix. I know of the following paper:

There, the author works on an algorithm of computing a adjugate matrix $\text{adj}(A)$ when $A$ is nearly singular or singular. For such matrices, one can make use of the factorizations. Suppose, we have found:

$$ A=XDY \tag{1} \label{eq1} $$ where $X$ and $Y$ are well conditioned, and $D$ is a diagonal matrix. Now, we can write the adjugate matrix, as follows:

$$ \text{adj}(A)=\text{det}(X)\text{det}(D)\text{det}(Y)\left(Y^{-1}D^{-1}X^{-1}\right) \label{eq2} \tag{2} $$

There are several standard decompositions that satisfy (with various guarantees) $\eqref{eq1}$: SVD, LU with full pivoting, pivoted QR and pivoted QLP. Now, the matrix $D$ enters the $\eqref{eq2}$ twice: as a $D^{-1}$ and $\text{det}(D)$ which seems like a problem in case the matrix is truly singular. The author of the paper argues (and justifies by the perturbation series analysis) that

  1. In floating-point arithmetic a true zero is unlikely
  2. If it really happens, a small perturbation should be applied to that and algorithm proceeds as per $\eqref{eq2}$ with perturbed zero entries of $D$.

The perturbation theory is unusual because although $\text{adj}(A)$ and $A^{-1}$ differ only by a scalar factor, the matrix $A^{-1}$ has singularities while $\text{adj}(A)$ is analytic - in fact, it is a multinomial in the elements of $A$. It turns out that multiplying by the determinant smooths out the singularities to give an elegant perturbation expansion.

...

However, if $A$ is ill-conditioned - that is, if $A$ is nearly singular - the inverse will be inaccurately computed. Nonetheless, we will show that this method, properly implemented, can give an accurate adjugate, even when the inverse has been computed inaccurately.

Take a look at the detailed discussion in the paper on the advantages and disadvantages of the proposed factorizations.

  • Going via SVD route guarantees the well-conditioning of $X$ and $Y$; however, finding their determinants might be tricky (even though, they are just signs: $\text{det}(X,Y)=\pm 1$.
  • On the contrary, both full-pivoted LU and pivoted QR should lead to an easy $\mathcal O(N^3)$ algorithm. For example full-pivoted LU:

$$ A=\Pi_\text{R} LDU\Pi_\text{C} $$

results in

$$ \text{adj}(A)=\text{det}(\Pi_\text{R})\text{det}(D)\text{det}(\Pi_\text{C})\left( \Pi_\text{C}^{T}U^{-1}D^{-1}L^{-1}\Pi_\text{R}^{T}\right) $$ where $\text{det}(\Pi_\text{R})=(-1)^{\text{number of row interchanges}}$ and all computations are straightforward.

So, that gives an $O(N^3)$ algorithm to compute the adjugate matrix since all the components are at most $O(N^3)$: finding the inverse of well-conditioned matrices, LU-decomposition, matrix-matrix multiplication, calculation of easy determinants. However, as opposed to SVD, the $X$ and $Y$ tend to be well-conditioned, but might not bee (see the detailed discussion in the paper). In practice, I don't think it would be an issue. And worst comes to worst, you might just have to use both methods in such special cases.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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