# Fast algorithm for computing cofactor matrix

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?

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)."?

• 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)? – GoHokies Jul 9 at 11:58
• 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. – ZUN LI Jul 9 at 14:36
• 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. – hardmath Jul 10 at 14:16

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.