Julia: How can we compute the adjoint or classical adjoint (Linear Algebra)?

Question

I wan to compute the classical adjoint in Julia 1.0

For this I copied the matrix given as an example in wikipedia

julia> B = [-3 2 -5; -1 0 -2; 3 -4 1]
3×3 Array{Int64,2}:
 -3   2  -5
 -1   0  -2
  3  -4   1

That seems to me to compute the transpose of B and not its adjoint. Instead, we should get this (from wikipedia):

and tried to get its adjoint using the adjoint() function which is mentionned in the Julia documentation here despite that the doc does not specifically said what this function does

julia> adjoint(B)
3×3 Adjoint{Int64,Array{Int64,2}}:
 -3  -1   3
  2   0  -4
 -5  -2   1

Instead I want to get this:

In Matlab I get indeed:

>> adjoint(B)

ans =

   -8.0000   18.0000   -4.0000
   -5.0000   12.0000   -1.0000
    4.0000   -6.0000    2.0000

Answer 1

Julia's adjoint is defined as the transpose of the complex conjugate of the input matrix. However, you seem to want the adjugate matrix:

The adjugate has sometimes been called the "adjoint", but today the "adjoint" of a matrix normally refers to its corresponding adjoint operator, which is its conjugate transpose.

You can compute the adjugate matrix by inverting, and then multiplying by the determinant:

julia> det(B) * inv(B)
3×3 Array{Float64,2}:
 -8.0  18.0  -4.0
 -5.0  12.0  -1.0
  4.0  -6.0   2.0

Thanks to @Antoine Levitt and @Syx Pek on the Julia Slack for giving the suggestion of inverting and multiplying by determinant.

---

Original answer:

The adjugate matrix seems to be the transpose of the matrix of cofactors. Below is a naïve implementation of finding cofactors:

# import Pkg; Pkg.add("InvertedIndices")
using InvertedIndices # for cleaner code, you can remove this if you really want to.
function cofactor(A::AbstractMatrix, T = Float64)
           ax = axes(A)
           out = similar(A, T, ax)
           for col in ax[1]
               for row in ax[2]
                   out[col, row] = (-1)^(col + row) * det(A[Not(col), Not(row)])
               end
           end
           return out
       end

Then, to find the adjugate, you only need to transpose (transpose(cofactor(B))).

The answer is:

julia> cofactor(B, Float64) |> transpose
3×3 Transpose{Float64,Array{Float64,2}}:
 -8.0  18.0  -4.0
 -5.0  12.0  -1.0
  4.0  -6.0   2.0

which is equivalent to what Matlab gives.

Edit: @Antoine Levitt on the Julia slack pointed out that this is essentially a rescaled inverse matrix, so if you figure out the scaling factor, you can just do inv(B) * scaling_factor (in the case of this matrix, it's 6).

Answer 2

Explicitly using the cofactor formula, as in the answer above, is very expensive — it has $O(n^5)$ complexity for an $n \times n$ matrix, because you are computing $n^2$ determinants of $(n-1) \times (n-1)$ matrices.

A more efficient $O(n^3)$ method is to use the formula $\mathrm{cofactor}(A) = \mathrm{adjugate}(A)^T = A^{-T} \det A$, but this fails if the matrix is singular (and may be very inaccurate for near-singular matrices). Instead, using the SVD $A = U \Sigma V^T$, one obtains $\mathrm{adjugate}(A) = V \mathrm{adjugate}(\Sigma) U^T \det(V^T U)$, and the adjugate of the diagonal matrix $\Sigma$ is easily computed.

A Julia implementation of the SVD-based $O(n^3)$ algorithm can be found in this forum post.