# Computing the adjugate/cofactors of a matrix

**URL:** <https://discourse.julialang.org/t/computing-the-adjugate-cofactors-of-a-matrix/125184>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [January 25, 2025, 12:06am UTC](https://discourse.julialang.org/t/computing-the-adjugate-cofactors-of-a-matrix/125184 "2025-01-25T00:06:04Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [January 25, 2025, 12:06am UTC](https://discourse.julialang.org/t/computing-the-adjugate-cofactors-of-a-matrix/125184/1 "2025-01-25T00:06:04Z")

</div>

For pedagogy (in our [matrix calculus course at MIT](https://github.com/mitmath/matrixcalc/)), I wanted to have a routine to compute the [adjugate of a matrix](https://en.wikipedia.org/wiki/Adjugate_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](https://scicomp.stackexchange.com/questions/33028/fast-algorithm-for-computing-cofactor-matrix).

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**
>
> ```julia
> 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](https://scicomp.stackexchange.com/questions/40345/calculate-determinant-of-unitary-matrices-based-on-svd-implementation), 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`.

---

<div class="post-metadata">

**Author:** ![Mervyn\_Larrier](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mervyn_larrier/32/37709_2.png) [@Mervyn\_Larrier](https://discourse.julialang.org/u/Mervyn_Larrier)\
**Post date:** [June 3, 2025, 12:04am UTC](https://discourse.julialang.org/t/computing-the-adjugate-cofactors-of-a-matrix/125184/2 "2025-06-03T00:04:19Z")

</div>

You should see if you can add it to the LinearAlgebra library as a fork, if possible. This works very well!