Calculate the inverse of a matrix

Question

Somehow, we don't yet have a challenge for finding the inverse of an arbitrarily-sized square matrix, despite having ones for 3x3 and 4x4, as well as a more complex version.

Your task is, given a square \$n\times n\$ non-singular matrix \$M\$, output the matrix \$M^{-1}\$ that satisfies

$$MM^{-1} = I_n$$

There are a number of methods and formulae for calculating \$M^{-1}\$, but one of the most well known is

$$M^{-1} = \frac1{\det(M)}\text{ adj}(M)$$

where \$\det\$ represents the determinant and \$\newcommand{\adj}{\text{adj}}\adj\$ the adjugate

Some definitions:

• \$I_n\$: The \$n\times n\$ identity matrix i.e. an \$n\times n\$ matrix where the leading diagonal consists entirely of \$1\$s and the rest \$0\$s
• Non-singular: the determinant of \$M\$ is guaranteed to be non-zero
• Determinant: a specific number that can be calculated for any given square matrix. Exact methods can be found in the Wikipedia article
• Adjugate: Formally, the transpose of the cofactor matrix of \$M\$. Informally, this is a operation on \$M\$ which takes determinants of submatrices in a specific way to construct a related matrix. Again, exact details can be found in the linked article.

For sake of simplicity, you may assume:

• The elements of \$M\$ will all be integers within the native bounds of your language
• \$n\$, nor \$n^2\$, will never exceed the maximum value in your language, and will always be greater than or equal to \$1\$
• The elements of \$M^{-1}\$ will never exceed the maximum value in your language (or minimum for negative values)
• \$M\$ will never be singular

No builtins are banned and you may use whatever (valid) method you like for calculating \$M^{-1}\$. It is acceptable if your program fails for some inputs due to floating point issues, so long as the underlying algorithm or method works for arbitrary matrices.

This is, of course, entirely optional, but if your answer consists entirely of a builtin, consider including a non-builtin method, simply for the sake of general interest.

Standard code-golf rules apply. This means you may input or output in any convenient format, and that standard loopholes are forbidden. The shortest code in bytes wins.

This script will take an input \$n\$ and generate an \$n\times n\$ matrix with random integers between \$-10\$ and \$10\$, along with it's inverse. You can use this for test cases.

---

Worked example

Lets take the \$3\times3\$ matrix \$M\$ as:

$$M = \left[\begin{matrix} 4 & -3 & 0 \\ -4 & -7 & 6 \\ 5 & 7 & 6 \end{matrix}\right]$$

We'll use the above formula, \$M^{-1} = \frac{\adj(M)}{\det(M)}\$ for this example.

First, we'll calculate \$\det(M)\$ by expanding along the third column:

$$\begin{align} \det(M) & = \left|\begin{matrix} 4 & -3 & 0 \\ -4 & -7 & 6 \\ 5 & 7 & 6 \end{matrix}\right| \\ & = 0\left|\begin{matrix} -4 & -7 \\ 5 & 7 \end{matrix}\right| - 6\left|\begin{matrix} 4 & -3 \\ 5 & 7 \end{matrix}\right| + 6\left|\begin{matrix} 4 & -3 \\ -4 & -7 \end{matrix}\right| \\ & = 0 - 6(4\cdot7 - -3\cdot5) + 6(4\cdot-7 - -3\cdot-4) \\ & = -6(28 + 15) + 6(-28 - 12) \\ & = -6\cdot43 + 6\cdot-40 \\ & = -498 \\ \therefore det(M) & = -498 \end{align}$$

We then need to calculate \$\adj(M)\$. As \$\adj(\cdot)\$ of a matrix is the transpose of the cofactor matrix, this essentially boils down to calculating the cofactor matrix of \$M\$, \$C_M\$:

$$\begin{align} \adj(M) & = C_M^T \\ & = \left[\begin{matrix} \left|\begin{matrix} -7 & 6 \\ 7 & 6 \end{matrix}\right| & \left|\begin{matrix} -4 & 6 \\ 5 & 6 \end{matrix}\right| & \left|\begin{matrix} -4 & -7 \\ 5 & 7 \end{matrix}\right| \\ \left|\begin{matrix} -3 & 0 \\ 7 & 6 \end{matrix}\right| & \left|\begin{matrix} 4 & 0 \\ 5 & 6 \end{matrix}\right| & \left|\begin{matrix} 4 & -3 \\ 5 & 7 \end{matrix}\right| \\ \left|\begin{matrix} -3 & 0 \\ -7 & 6 \end{matrix}\right| & \left|\begin{matrix} 4 & 0 \\ -4 & 6 \end{matrix}\right| & \left|\begin{matrix} 4 & -3 \\ -4 & -7 \end{matrix}\right| \end{matrix}\right]^T \\ & = \left[\begin{matrix} -84 & 54 & 7 \\ 18 & 24 & -43 \\ -18 & -24 & -40 \end{matrix}\right]^T \\ & =\left[\begin{matrix} -84 & 18 & -18 \\ 54 & 24 & -24 \\ 7 & -43 & -40 \end{matrix}\right] \end{align}$$

Finally, having calculated both \$\det(M)\$ and \$\adj(M)\$, we divide each element of \$\adj(M)\$ by \$\det(M)\$ to compute the final output, \$M^{-1}\$:

$$\begin{align} M^{-1} & = \frac{\adj(M)}{\det(M)} \\ & = \left[\begin{matrix} \frac{-84}{-498} & \frac{ 18}{-498} & \frac{-18}{-498} \\ \frac{ 54}{-498} & \frac{ 24}{-498} & \frac{-24}{-498} \\ \frac{ 7}{-498} & \frac{-43}{-498} & \frac{-40}{-498} \end{matrix}\right] \\ & = \left[\begin{matrix} \frac{ 14}{ 83} & \frac{-3}{ 83} & \frac{ 3}{ 83} \\ \frac{ -9}{ 83} & \frac{-4}{ 83} & \frac{ 4}{ 83} \\ \frac{ -7}{498} & \frac{43}{498} & \frac{20}{249} \end{matrix}\right] \end{align}$$

Alternatively, as decimals, \$M^{-1}\$ is

[[ 0.1686746987951807,   -0.03614457831325301, 0.03614457831325303],
 [-0.10843373493975902,  -0.04819277108433735, 0.04819277108433734]
 [-0.014056224899598388,  0.08634538152610442, 0.08032128514056225]]

Answer 1

Octave, 57 bytes

A=input('');V=A'/trace(A*A');for i=1:1e4V=2*V-V*A*V;end
V

Try it online!

This is not particularly well golfed, but I wanted to advertise an approach that could be useful for other non-builtin answers.

This uses the Hotelling-Bodewig scheme:

$$ V_{i+1} = V_i\left(2I - AV_i\right)$$

Which iteratively computes the inverse of a non singular matrix. This is guaranteed to converge for \$\left\lVert I - AV_0\right\rVert < 1\$ (under a suitable matrix norm). Choosing the \$V_0\$ is difficult, but Soleymani, F. shows in "A New Method For Solving Ill-Conditioned Linear Systems" that the inital guess \$V_0 = \frac{A^T}{\text{tr}(AA^T)}\$ will always satisfy this condition, so the system is numerically stable.

What makes this a particularly attractive approach to other potential answers is that we don't require any builtin determinant or inverse functions. The most complex part is just matrix multiplication, since the transpose and trace are trivial to compute.

I have chosen 1e4 iterations here to make the runtime somewhat reasonable, although you could of course push it to 1e9 with no loss of byte count.

---

-10 thanks to xnor for noting we don't need to construct an identity matrix.

Answer 2

APL (Dyalog Unicode), 1 byte^SBCS

⌹

Try it online!

The domino primitive is a very interesting APL "built-in". It already featured in another 1-byte answer of my own where it was used to solve a least-squares problem. When applied to a square matrix, ⌹ tries to find the matrix inverse of its argument.

Many golfing languages will also have a built-in for this... But mind you, APL is not a golfing language, although it is terse enough to be very competitive and, in cases like this, win.

Answer 3

R, 5 bytes

solve

Try it online!

Nothing new here... Basically, the code solve(A, B) solves \$AX = B\$, but when \$B\$ is not given, it is treated as identity matrix, thus giving us the inverse as the result.

Answer 4

R, 72 61 bytes

function(A,V=t(A/sum(diag(A%*%t(A))))){for(i in 1:1e4)V=2*V-V%*%A%*%V;V}

Try it online!

Porting Sisyphus' answer is not futile at all...and thanks to Sisyphus for -11 bytes.

Observes that \$Tr(AA^T)=\sum\limits_{i,j}a_{ij}^2\$.

R, 94 bytes

function(M)outer(k<-1:dim(M),k,Vectorize(function(j,i)det(M[-i,-j,drop=F])*(-1)^(i+j)))/det(M)

Try it online!

Thanks to Robin Ryder for fixing a bug and making this actually work.

Calculates \$A^{-1}\$ using the adjugate/determinant method.

Answer 5

Python 2, 228 bytes

from random import*
a=input()
exec"""$:j,J=i,I;J+=[j==i $]
while~-all(I[i]$):shuffle(a)
$:
 j,J=i,I
 $:
	if j-i:I[:]=[y-I[j]*x/J[j]for x,y in zip(J,I)]
$:print[x/I[i]for x in I][len(a):]""".replace("$","for i,I in enumerate(a)")

Try it online!

Augment the matrix with the identity matrix, then apply Gauss–Jordan elimination. I don't know if this is the shortest approach, but it's the one I wanted to try golfing down.

I use while not all(a[i][i]for i in r):shuffle(a) to move zeros off the diagonal. This loop will definitely terminate, because if there is no permutation of the rows of \$A\$ that makes the diagonal free of zeros, then \$\det(A)=0\$, which we are guaranteed is not the case. This can be seen from the Leibniz formula for \$\det(A)\$:

$$\det(A) = \sum_{\sigma \in S_n} \text{sgn}(\sigma) \prod_{i=1}^n a_{\sigma(i),i}$$

“There is no permutation \$\sigma\$ of the rows that makes the diagonal free of zeros” can be equivalently rephrased as “\$\prod_{i=1}^n a_{\sigma(i),i}\$ is always 0, for all \$\sigma\$” which causes this whole formula to be 0.

Answer 6

APL (Dyalog Unicode), 27 23 bytes

-4 bytes thanks to Adám!

This implements the method advertised by Sisyphus.

⊢(⊢+⊢-⊢+.×+.×)⍣≡⍉÷,+.×,

Try it online!

A function that takes the matrix as the right argument.

+.× computes the dot product (sum of element-wise product) of the flattened matrix , and the flattened matrix ,. This is sum of the squared matrix entries, which is equal to \$tr(AA^T)\$.
⍉÷ divides the transposed matrix by the trace.

(⊢+⊢-⊢+.×+.×)⍣≡ is a function which is applied to the original matrix \$A\$ as a left argument (⊢) and \$A^T\div tr(AA^T)\$ as the right argument (⍉÷,+.×,):

⊢+⊢-⊢+.×+.× takes the current matrix \$V\$ on its right and the input matrix \$A\$ on its left and executes one iteration step:
+.× is the inner product of + and ×. Given two matrices, this calculates their product. In this case \$ A \times V \$.
⊢ is the right argument \$V\$, ⊢+.× the product \$V \times (A \times V)\$.
⊢- subtracts this from the right argument: \$V-V \times A \times V\$.
⊢+ adds this to the right argument: \$V+V-V \times A \times V\$.

⍣≡ applies the function on its left until the result doesn't change. Because of the way equality testing works in Dyalog APL, this actually terminates.

Answer 7

Julia 1.0, 3 bytes

inv

Try it online!

Yet another short built-in solution.

Answer 8

JavaScript (ES6), 169 bytes

This computes \$M^{-1} = \dfrac{\operatorname{adj}(M)}{\det(M)}\$

M=>M.map((r,y)=>r.map((_,x)=>D(h(M,x).map(r=>h(r,y)))*(x+y&1?-1:1)/D(M)),h=(a,n)=>a.filter(_=>n--),D=M=>+M||M.reduce((s,[v],i)=>s+(i&1?-v:v)*D(h(M,i).map(r=>h(r,0))),0))

Try it online!

Answer 9

J, 2 bytes

%.

Try it online!

Same as APL, but more powerful, as J can produce exact rational matrix when given a matrix of extended integers as input.

Answer 10

05AB1E, 38 22 21 20 bytes

˜nO/øтF©I2Føδ*O®}·s-

Port of @Sisyphus' Octave answer, so make sure to upvote him!!
-16 bytes thanks to @ovs.

Try this online.

Code explanation:

˜            # Flatten the (implicit) input-matrix to a single list
 n           # Square each value in this list
  O          # Take the sum (this is the trace of M*M')
   /         # Divide each value in the (implicit) input-matrix by this trace
    ø        # Zip/transpose this matrix; swapping rows/columns
тF           # Loop 100 times:
  ©          #  Store the current matrix in variable `®` (without popping)
   I         #  Push the input-matrix
    2F       #  Loop 2 times:
      ø      #   Zip/transpose the top matrix; swapping rows/columns
       δ     #   Apply double-vectorized with the top two matrices:
        *    #    Multiply
         O   #   Sum each inner row
          ®  #   Push the matrix from variable `®` again
     }·      #  After the inner loop: double all values in matrix `®`
       s     #  Swap so the calculated matrix VMV is at the top again
        -    #  Subtract this VMV from the 2V
             # (after the outer loop, the resulting matrix is output implicitly) 

---

Original answer (38 bytes) and detailed explanation:

εUεX*O]Å\OIøs/тFxs©εUIøεX*O}U®øεX*O}}-

Try it online.

05AB1E has barely any useful builtins for matrices, not even matrix multiplication. So almost everything has to be done manually..

Since I'm an absolute noob in math, I'm gonna explain everything in full detail to help others like me who want to do this challenge without any builtins, and also to keep this answer self-contained.

Step 1) Matrix multiplication of the input-matrix \$M\$ with it's transpose: \$M\times M'\$:

If we have a matrix \$A\$ and \$B\$ and want to do matrix-multiplication \$AB\$, we take the dot-product of every \$i^{th}\$ row of \$A\$ and \$j^{th}\$ column of B for every coordinate \$i,j\$ in the two matrices.

For example, if we use the matrix in the challenge description:

\$M = \left[\begin{matrix} 4 & -3 & 0 \\ -4 & -7 & 6 \\ 5 & 7 & 6 \end{matrix}\right]\$

We can for example calculate the values in the top row of the resulting \$M\times M'\$ matrix with:

Top-left: \$4\times4+-3\times-3+0\times0 = 25\$
Top-center: \$4\times-4+-3\times-7+0\times6=5\$
Top-right: \$4\times5+-3\times7+0\times6 = -1\$

I've done matrix multiplication in 05AB1E before in this answer of mine, so I've used that code snippet here as well. Since we want to multiply the input-matrix by it's transpose, we actually won't need the transpose builtin here.

ε                  # Map over each row of the (implicit) input-matrix
 U                 #  Pop and store the current row in variable `X`
  ε                #  Map over each row of the (implicit) input-matrix again
   X*              #   Multiply the values of the current row by the values at the same
                   #   positions in row `X`
     O             #   And take the sum of this row
]                  # Close both maps

Try just this step online.

Step 2) Take the trace of this new matrix: \$(M\times M')^T\$

The trace of a square matrix is basically the sum of its main diagonal (the values of the top-left to the bottom-right).

Å\                 # Take the main diagonal of the matrix of step 1
  O                # And sum the values in this list together

Try the first two steps online.

Step 3) Divide all values in the transposed matrix by this trace we calculated:

I                  # Push the input-matrix
 ø                 # Zip/transpose it; swapping rows/columns
  s                # Swap so the trace we calculated it at the top of the stack
   /               # And divide each value in the transposed matrix by this trace

Try the first three steps online.

Step 4) Repeat the following steps (5 through 8) enough times for the answer to not change anymore:

Since this program isn't very fast in 05AB1E, I've decided to loop just 100 times, but this can be increased to improve the accuracy of the decimal results (I've verified with @Sisyphus' Octave answer that changing the 1e4 to 1e2 still holds the same result for most matrices).

тF                 # Loop 100 times:

I'm not sure if the values will eventually not change anymore if we loop enough times. If this is the case we could (in theory) save a byte by changing this тF to Δ (loop until the result no longer changes).

(Let's call the intermediate matrix inside this loop \$V\$ for the explanations of the following steps.)

Step 5) Double each value in the current matrix: \$2V\$:

  x                #  Double each value in the current matrix V (without popping)

Try the first five steps online, excluding the loop of step 4.

Step 6) Do matrix multiplication again for \$VM\$ (where \$M\$ is the input-matrix):

   s               #  Swap to take the non-doubled matrix V at the top again
    ©              #  Store this matrix V in variable `®` (without popping)
     ε             #  Map over each row of matrix V:
      U            #   Pop the current row, and store it in variable `X`
       I           #   Push the input-matrix M
        ø          #   Zip/transpose; swapping rows/columns
         ε         #   Map over each row of this transposed matrix M':
          X*       #    Multiply the values in the current row by row `X`
            O      #    And take the sum

Try the first six steps online, excluding the loop of step 4.

Step 7) And do matrix multiplication yet again right after: \$VMV\$:

         }         #   Close the inner map
          U        #   Pop and store this as new `X`
           ®       #   Push the matrix V from variable `®`
            ø      #   Zip/transpose; swapping rows/columns
             ε     #   Map over each row of this transposed matrix V':
              X*   #    Multiply the values in the current row by row `X`
                O  #    And take the sum
     }}            #  Close both the inner and outer maps

Try the first seven steps online, excluding the loop of step 4.

Step 8) Subtract the values at the same positions of these two matrices from one another: \$2V-VMV\$:

       -           #  Subtract matrix VMV from 2V

Try the first eight steps online, excluding the loop of step 4.

And after the loop is done, the resulting matrix is output implicitly.

Answer 11

Python 2, 188 bytes

lambda a:[[c(a,j,i)/d(a)for j,_ in e(a)]for i,_ in e(a)]
c=lambda a,i,j:(-1)**(i+j)*d([b[:j]+b[j+1:]for I,b in e(a)if i-I])
d=lambda a:a==[]or sum(b[0]*c(a,i,0)for i,b in e(a))
e=enumerate

Try it online!

The top lambda computes \$A^{-1} = \frac{1}{\det(A)}\text{adj}(A)\$.

d(a) computes the determinant and c(a,i,j) computes cofactors.

Answer 12

Excel, 29 bytes

=MINVERSE(OFFSET(A2,,,A1,A1))

Straightforward application of the MINVERSE() function. It's boring but I got excited about Excel having a built-in for something. Input \$n\$ in A1, the matrix starting in A2, and the formula anywhere the spill won't interfere.

Answer 13

Scala, 237 232 bytes

Uses the method from Sisyphus's answer. Go upvote that!

m=>{val h=m.indices
Seq.iterate(m.transpose.map(_.map(_/m.flatten.map(x=>x*x).sum)),9999){v=>h.map(i=>h.map{j=>2*v(i)(j)-(h.map(k=>v(i).zip(m.transpose.apply(k))map(t=>t._1*t._2)sum),v.transpose.apply(j)).zipped.map(_*_).sum})}last}

Try it online!

h is just a range from 0 until n to reuse later (mostly because Scala doesn't have matrix multiplication builtins). The function makes a sequence of 9999 elements and takes the last element. The first element is the transpose of m divided by the trace of m times its transpose. Subsequent elements are calculated with 2*v-v*m*v, where v was the previous element.

To calculate \$V_0\$ (It turns out the trace of m times its transpose is just the sum of squares of all of m's cells):

m.transpose.map(            //For every row in m's transpose
  _.map(                    //For every cell in that row
    _ /                     //Divide it by (trace(M * M's transpose))
      m.flatten             //Turn m into a 1D list
        .map(x=>x*x)        //Square each cell
        .sum))              //Add them up

To calculate subsequent elements, we use \$2V - (VA)V\$, but you have to map over h instead of over v itself:

h.map(i =>                 //For every i in [0, n)
  h.map{j =>               //For every j in [0, n)
    2*v(i)(j) -            //2V at these coordinates minus
    <(v * m * v)[i][j]> }) //v*m*v at these coordinates (see explanation below)

To calculate (v*m)[i]:

h.map(k =>                //k is the index of a row in [0, n)
  v(i).zip(               //Zip column i of v with
    m.transpose.apply(k)  //Row k of m (apply is used for indexing here)
  ) map(t=>t._1*t._2)     //Multiply v(i)(j) with m(k)(i)
  sum                     //Add then up
)

And getting the cross product of that with row j of v uses pretty much the same approach.

---

Scala, 346 342 bytes

Saved 4 bytes thanks to @corvus_192!

type M=Seq[Seq[Double]]
def c(m:M)={val I=m.indices;I.map(i=>I.map(j=>m(i)(j)*math.pow(-1,i+j)))}
def d(m:M):(M,Double)=if(m.size<2)m->m(0)(0)else{val I=m.indices
val M=I.map(i=>I.map{j=>d(I.filter(i!=_)map(k=>I.filter(j!=_)map(m(k))))._2})
c(M)->c(m).head.zip(M.head).map(t=>t._1*t._2).sum}
def i(m:M)=d(m)._1.transpose.map(_.map(_/d(m)._2))

Try it in Scastie!

As you can see, I'm not very good at math.

Answer 14

Wolfram Language (Mathematica), 7 bytes

Inverse

Try it online!

Answer 15

Jelly, 3 bytes

æ*-

Try it online.

Explanation:

     # Full program taking a single integer-matrix as argument
æ*   #  Matrix exponentiation
  -  #  with -1
     #  (after which the result is output implicitly)

Answer 16

Matlab ⁶ 3 bytes

inv

Computes and prints the inverse of a square matrix. Pretty boring in-built solution. Thanks to @Bubbler for the clarification and -3 bytes.

Answer 17

Charcoal, 48 bytes

≔Ｅθ∕Ｅθ§λκΣＥＸθ²ΣληＦφＵＭηＥκ⁻⊗μΣＥθ×ΣＥθקκς§ρπ§§ηπνＩη

Try it online! Link is to verbose version of code. Explanation: Another port of @Sisyphus's answer.

≔Ｅθ∕Ｅθ§λκΣＥＸθ²Σλη

Transpose the input and divide it by the sum of squares of all the elements. Sadly neither sum nor divide fully vectorise, so I have to divide a row at a time and calculate the sum via a nested loop.

Ｆφ

Repeat 1000 times, which should be enough for floating-point precision.

ＵＭηＥκ⁻⊗μΣＥθ×ΣＥθקκς§ρπ§§ηπν

Calculate the matrix multiplication and subtraction in-place. Charcoal doesn't have any vector or matrix operations, so we have to loop over the rows and columns manually, but there are a couple of places where we can share variables which saves us a couple of bytes each.

Ｉη

Output the array. (Note that each element is output on its own line and each row is double-spaced from the previous.)

Answer 18

MATL, 4 bytes

-1Y^

Try it online!

Explanation

-1Y^
-1   : Push -1 onto the stack
  Y^ : Raise implicit input to -1 power

Answer 19

PARI/GP 8 bytes

f(A)=1/A

Built-in function; particularly efficient for matrices with integer matrix elements.

Most likely, the method used is similar to the one described in this refererence: H. Haramotu, M. Matsumotu, A p-adic algorithm for computing the inverse of integer matrices, Journal of Computational and Applied Mathematics, Volume 225, Issue 1, 1 March 2009, Pages 320-322.

Answer 20

SageMath, ^{14 13} 11 bytes

Saved a byte thanks to FryAmTheEggman!!!
Saved 2 bytes thanks to Sisyphus!!!

lambda M:~M

Try it online!

Inputs any square matrix and returns its inverse.

Answer 21

Ruby -rmatrix, ²³ 19 bytes

->a{Matrix[*a].inv}

Try it online!

Returns the result as a Ruby matrix object.

-4 bytes from Dingus.

Answer 22

q, 3 bytes

inv

k, 1 byte

!

Does exactly as requested. Input must be a float matrix

Answer 23

Fortran (GFortran), 350 376 bytes

real,allocatable,dimension(:,:)::A,B,U;read*,n;allocate(A(n,n));B=A;U=A;read*,((A(i,j),j=1,n),i=1,n)
do i=1,n;do j=1,n;U(i,j)=merge(1,0,i==j);B(i,j)=merge(1.,1/A(i,j),A(i,j)==0);enddo;enddo
do 4 l=1,99;do 4 i=1,n;do 4 j=1,n;s=U(i,j);do k=1,n;m=i;if(k/=m)s=s-A(i,k)*B(k,j);enddo;B(m,j)=s/A(i,m)
4 continue;print'(3(x,f9.6))',((B(i,j),j=1,n),i=1,n)
end

Try it Online!   _{376 bytes}

Adapted from here. Saved 26 bytes using merge instead of if.