Computing the adjugate/cofactors of a matrix

Specific Domains Numerics
1

For pedagogy (in our matrix calculus course at MIT), I wanted to have a routine to compute the adjugate of a matrix, since this arises in the gradient of the determinant. This is defined by:

\mathrm{adj}(A) = A^{-1} \det A = \mathrm{cofactor}(A)^T

Computing it using det(A) * inv(A) is problematic, however, because it fails if the matrix is singular (and could be numerically unstable if the matrix is badly conditioned).

Instead, one can compute it using an SVD.

For reference, I’ve included the code below. It might be worth putting in a package, although the adjugate does not seem very practical by itself (it quickly over/underflows for large matrices, similar to the determinant).

Adjugate implementation 
using LinearAlgebra

function adjugate(A::AbstractMatrix)
    ishermitian(A) && return adjugate(Hermitian(A))
    LinearAlgebra.checksquare(A)
    F = svd(A)
    return F.V * (adjugate(Diagonal(F.S)) * det(F.Vt * F.U)) * F.U'
end

function adjugate(A::LinearAlgebra.RealHermSymComplexHerm)
    F = eigen(A)
    return F.vectors * adjugate(Diagonal(F.values)) * F.vectors'
end

function adjugate(D::Diagonal{<:Number})
    d = D.diag
    Base.require_one_based_indexing(d)
    length(d) < 2 && return Diagonal(one.(d))

    # compute diagonal of adj(D): dadj[i] = prod(d) / d[i], but avoid dividing by zero
    dadj = similar(d)
    prod = one(eltype(d))
    for i = 1:length(d)
        dadj[i] = prod
        prod *= d[i]
    end
    prod = one(eltype(d))
    for i = length(d):-1:1
        dadj[i] *= prod
        prod *= d[i]
    end
    
    return Diagonal(dadj)
end

PS. Part of the SVD-based adjugate algorithm requires one to compute the determinant of the U and V factors. Annoyingly, even though in principle these determinants are implicit in the SVD algorithm, there appears to be no efficient way to access that information, forcing one to just call det. On the other hand, the det computation is probably negligible compared to svd.

5 Likes