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

Improve this question
$\endgroup$
3
  • 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, 2019 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, 2019 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, 2019 at 14:16

1 Answer 1

Reset to default
7
$\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 an 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 the 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 $\mathcal O(N^3)$ algorithm to compute the adjugate matrix since all the components are at most $\mathcal 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 be (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.

Improve this answer
$\endgroup$
3
  • $\begingroup$ Since $\det(D) = \prod_{i=1}^n d_{ii}$, it seems to me that one can use the identity $\det(D) D^{-1} = \operatorname{diag}(e_i)$, where $e_i = \prod_{j\neq i} d_{jj}$. This gives an algorithm that works without divisions by zero even if $D$ is exactly singular: no need to introduce "small perturbations". $\endgroup$
    – Federico Poloni
    Feb 23, 2020 at 20:09
  • 1
    $\begingroup$ @FedericoPoloni interesting, I did not think of it while reading that paper. That seems like a small but valuable improvement over the 1998 paper to me. $\endgroup$
    – Anton Menshov ♦
    Feb 23, 2020 at 20:57
  • $\begingroup$ Actually that formula is literally in the Stewart paper (Eqns 1.4 and 1.5) for the special case of the SVD, but surprisingly the author does not use it later and resorts to that perturbation argument. Maybe I am missing something here, I did not read the paper carefully. $\endgroup$
    – Federico Poloni
    Feb 23, 2020 at 21:31

Your Answer

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

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