Find the Matrix Power

Question

Problem

Create a program or function that can calculate the result of a matrix raised to the n^th power. Your code will take an arbitrary square matrix A and a non-negative integer n, and return a matrix with the value Aⁿ.

Restrictions

Built-in functions that compute the matrix power and matrix product are not allowed.

The rest of the standard rules for code-golf apply.

Explanation

Given a square matrix A, the value of Aⁿ = A A ⋯ A (repeated matrix product of A with itself, n times). If n is positive, the standard just mentioned is used. When n is zero, the identity matrix with the same order of A is the result.

Goal

This is code-golf and the shortest code wins.

Test Cases

Here, A is the input matrix, n is the input integer, and r is the output matrix where r = Aⁿ.

n = 0
A = 62 72
 10 34
r = 1 0
 0 1
n = 1
A = 23 61 47
 81 11 60
 42 9 0
r = 23 61 47
 81 11 60
 42 9 0
n = 2
A = 12 28 -26 3
 -3 -10 -13 0
 25 41 3 -4
 -20 -14 -4 29
r = -650 -1052 -766 227
 -331 -517 169 43
 332 469 -1158 -53
 -878 -990 574 797
n = 4
A = -42 -19 18 -38
 -33 26 -13 31
 -43 25 -48 28
 34 -26 19 -48
r = -14606833 3168904 -6745178 4491946
 1559282 3713073 -4130758 7251116
 8097114 5970846 -5242241 12543582
 -5844034 -4491274 4274336 -9196467
n = 5
A = 7 0 -3 8 -5 6 -6
 6 7 1 2 6 -3 2
 7 8 0 0 -8 5 2
 3 0 1 2 4 -3 4
 2 4 -1 -7 -4 -1 -8
 -3 8 -9 -2 7 -4 -8
 -4 -5 -1 0 5 5 -1
r = 39557 24398 -75256 131769 50575 14153 -7324
 182127 19109 3586 115176 -23305 9493 -44754
 146840 31906 -23476 190418 -38946 65494 26468
 42010 -21876 41060 -13950 -55148 19290 -406
 44130 34244 -35944 34272 22917 -39987 -54864
 1111 40810 -92324 35831 215711 -117849 -75038
 -70219 8803 -61496 6116 45247 50166 2109

• Dennis: What about built-ins for matrix product or matrix inversion?
• miles: @Dennis I was considering banning those also, but it felt too restrictive.
• xnor: For languages without built-in matrix inversion, this strike me as a chameleon challenge because inverting a matrix from scratch seems much harder than taking the iterated product.
• Luis Mendo: I agree with @xnor. And what if a language doesn't have matrix inversion but has matrix power? Can A^-1 be used as a substitute for inv(A)?
• xnor: Is exp(k*log(M)) allowed? (Though it might not work because of non-unique branches.)

Answer 1

Jelly, (削除) 17 (削除ここまで) (削除) 16 (削除ここまで) 15 bytes

×ばつ'3S€€
LṬ€z0Ç¡

Try it online!

Permalinks with grid output: test case 1 | test case 2 | test case 3 | test case 4 | test case 5

How it works

LṬ€z0Ç¡ Main link. Arguments: A (matrix), n (power)
L Get the length l of A.
 Ṭ€ Turn each k in [1, ..., l] into a Boolean list [0, 0, ..., 1] of length k.
 z0 Zip; transpose the resulting 2D list, padding with 0 for rectangularity.
 This constructs the identity matrix of dimensions k×ばつk.
 Ç¡ Execute the helper link n times×ばつ'3S€€ Helper link. Argument: B (matrix)
Z Zip; transpose rows and columns of B.
 3 Yield A.
 ×ばつ' Spawned multiplication; element-wise mutiply each rows of zipped B (i.e.,
 each column of B) with each row of A.
 This is "backwards", but multiplication of powers of A is commutative.
 S€€ Compute the sum of each product of rows.

Answer 2

MATL, 20 bytes

XJZyXyi:"!J21ドル!*s1$e

Try it online!

Explanation

This avoids the matrix multiplication by doing it manually, using element-wise multiplication with broadcast followed by vectorized sum. Specifically, to multiply matrices M and N, both of size s×ばつs:

1. Transpose M. Call the resulting matrix P.
2. Permute the dimensions of N such that N is "turned" with a rotation axis along the first dimension, giving an s×ばつs 3D array, say Q.
3. Multiply each element of P times each element of Q, with implicit broadcast. This means that P is automatically replicated s times along the third dimension, and Q is replicated s times along the second, to make them both s×ばつs×ばつs, before the actual element-wise multiplication takes place.
4. Sum along the first dimension to yield a ×ばつs×ばつs array.
5. Squeeze the leading singleton dimension out, to produce an s×ばつs result.

Commented code:

XJ % take matrix A. Copy to clipboard J
ZyXy % generate identity matrix of the same size
i: % take exponent n. Generate range [1 2 ... n] (possibly empty)
" % for each element in that range
 ! % transpose matrix with product accumulated up to now (this is step 1)
 J % paste A
 21ドル! % permute dimensions: rotation along first-dimension axis (step 2)
 * % element-wise multiplication with broadcast (step 3)
 s % sum along first dimension (step 4)
 1$e % squeeze out singleton dimension, i.e. first dimension (step 5)
 % end for. Display

• CalculatorFeline: Fails for 0....
• Luis Mendo: @CatsAreFluffy Thanks! Corrected

Answer 3

APL, (削除) 32 (削除ここまで) 31 chars

{⍺=0:(⍴⍵)×ばつ⍨⍣(⍺-1)⊣⍵}

Left argument the power to raise to, right argument the matrix. The hardest (most space consuming) bit is building the identity matrix for the case where the desired exponent is 0. The actual computation is based on the fact that the generalised inner product (.) with + and ×ばつ as operands is effectively the matrix product. This combined with the power ⍣ operator ("repeat") forms the meat of the solution.

• Adalynn: 1: Are you THE Stefano that beat Dan&Nick by one byte in the 2016 year game?! 2. (1+≢⍵)↑1 => 1↑⍨1+≢⍵ to save one byte.
• lstefano: Yes, that's me.

Answer 4

K (ngn/k), 14 bytes

{y(+/x*)'/=#x}

Try it online!

-2 bytes thanks to @coltim on the k tree. The inner train multiplies x on the right side of the current matrix (instead of left side).

Why (+/x*)' is also a matmul:

(+/(e f;g h)*)' (a b;c d)
( (+/(e f;g h)*) a b; (+/(e f;g h)*) c d )
( +/ (ae af;bg bh) ; +/ (ce cf;dg dh) )
((ae+bg) (af+bh); (ce+dg) (cf+dh))

K (ngn/k), 16 bytes

{y(+/'x*\:)/=#x}

Try it online!

Original solution. Takes [matrix; n] and computes matrix^n.

How it works

{y(+/'x*\:)/=#x} x:matrix; y:n
 =#x Identity matrix of size of x
 y( )/ Repeat y times:
 x*\: Multiply each row of x to each column of current matrix
 +/' Sum each to get the matrix product

Why +/'*\: is a matmul:

(a b;c d) *\: (e f;g h)
multiply-eachleft expands to
(a b * (e f;g h) ; c d * (e f;g h))
atomic application gives
((ae af
 bg bh)
 (ce cf
 dg dh))
so applying +/' on this will give the result of matmul.

Answer 5

Sage, 112 bytes

lambda A,n:reduce(lambda A,B:[[sum(map(mul,zip(a,b)))for b in zip(*B)]for a in A],[A]*n,identity_matrix(len(A)))

Try it online

Explanation:

The inner lambda (lambda A,B:[[sum(map(mul,zip(a,b)))for b in zip(*B)]for a in A]) is a straightforward implementation of matrix multiplication. The outer lambda (lambda A,n:reduce(...,[A]*n,identity_matrix(len(A)))) uses reduce to compute the matrix power by iterated matrix multiplication (using the aforementioned homemade matrix multiplication function), with the identity matrix as the initial value to support n=0.

Answer 6

Julia, (削除) 90 (削除ここまで) (削除) 86 (削除ここまで) 68 bytes

f(A,n)=n<1?eye(A):[A[i,:][:]⋅f(A,n-1)[:,j]for i=m=1:size(A,1),j=m]

This is a recursive function that accepts a matrix (Array{Int,2}) and an integer and returns a matrix.

Ungolfed:

function f(A, n)
 if n < 1
 # Identity matrix with the same type and dimensions as the input
 eye(A)
 else
 # Compute the dot product of the ith row of A and the jth column
 # of f called on A with n decremented
 [dot(A[i,:][:], f(A, n-1)[:,j]) for i = (m = 1:size(A,1)), j = m]
 end
end

Try it online! (includes all but the last test case, which is too slow for the site)

Saved 18 bytes thanks to Dennis!

Answer 7

Python 2.7, (削除) 158 (削除ここまで) 145 bytes

The worst answer here, but my best golf in Python yet. At least it was fun learning how to do matrix multiplication.

Golfed:

def q(m,p):
 r=range(len(m))
 if p<1:return[[x==y for x in r]for y in r]
 o=q(m,p-1)
 return[[sum([m[c][y]*o[x][c]for c in r])for y in r]for x in r]

Explanation:

#accepts 2 arguments, matrix, and power to raise to
def power(matrix,exponent):
 #get the range object beforehand
 length=range(len(matrix))
 #if we are at 0, return the identity
 if exponent<1:
 #the Boolean statement works because Python supports multiplying ints by bools
 return [[x==y for x in length] for y in length]
 #otherwise, recur
 lower=power(matrix,exponent-1)
 #and return the product
 return [[sum([matrix[c][y]*lower[x][c] for c in length]) for y in length] for x in length]

Answer 8

Jelly, 11 bytes

L=þ`Zḋþ3ƊƓ¡

Try it online!

A more modern update to Dennis' answer, be sure to give that an upvote.

Additionally, this is a 9 byte answer that takes the dimensions of the matrix as the first 2 command line args and the matrix as the third.

Both take the power via STDIN.

How it works

L=þ`Zḋþ3ƊƓ¡ - Main link. Takes A on the left
L - Length of A, L
 ` - Using L as both arguments:
 þ - Create an LxL matrix, where the element at (i,j) is:
 = - Does i = j?
 This creates an identity matrix LxL
 Ɠ - Read an integer from STDIN, n
 Ɗ ¡ - Do the following n times:
 Z - Transpose
 3 - Yield A
 þ - Pair all rows of the transposed matrix with the rows of A and, over each row:
 ḋ - Dot product

• Bubbler: I was about to ask for this (an updated Jelly answer) when coltim found the 14-byter in K, but you ninjad me :P

Answer 9

Julia, (削除) 27 (削除ここまで) 24

a$n=round.(exp(log(a)n))

I am not sure what is allowed, but...

julia> a×ばつ5 Matrix{Int64}:
 35 18 40 37 77
 31 5 45 23 73
 62 67 29 85 97
 20 9 83 70 65
 2 13 53 59 52
julia> round.(exp(log(a)5)) ≈ a^5
true

24 thanks to @MarcMush.

Answer 10

JavaScript (ES6), 123 bytes

(n,a)=>[...Array(n)].map(_=>r=m(i=>m(j=>m(k=>s+=a[i][k]*r[k][j],s=0)&&s)),m=g=>a.map((_,n)=>g(n)),r=m(i=>m(j=>+!(i^j))))&&r

I had a 132 byte version using reduce but I was just mapping over a so often that it turned out to be 9 bytes shorter to write a helper function to do it for me. Works by creating the identity matrix and multiplying it by a longhand n times. There are a number of expressions that return 0 or 1 for i == j but they all seem to be 7 bytes long.

• user100690: (n,a)=>... -> n=>a=>...

Answer 11

Python 3, 147 bytes

def f(a,n):
 r=range(len(a));r=[[i==j for j in r]for i in r]
 for i in range(n):r=[[sum(map(int.__mul__,x,y))for y in zip(*a)]for x in r]
 return r

Try it online!

Answer 12

R, 49 bytes

f=function(m,n)`if`(n,m%*%f(m,n-1),diag(nrow(m)))

Recursive function that takes a matrix and the power n to raise it to. Recursively calls %*%, which computes the dot-product. The initial value for the recursion is the identity matrix of the same size as m. Since m %*% m = m %*% m %*% I, this works just fine.

Answer 13

Python 3, 128 bytes

def f(A,n):r=range(len(A));return n and[[sum(A[i][j]*x[i]for i in r)for j in r]for x in f(A,n-1)]or[[i==j for i in r]for j in r]

Try it online!

Answer 14

Python 2, 131 bytes

f=lambda M,n:n and[[sum(map(int.__mul__,r,c))for c in zip(*f(M,n-1))]for r in M]or[[0]*i+[1]+[0]*(len(M)+~i)for i in range(len(M))]

Try it online!

Answer 15

05AB1E, 11 bytes

2āDδQ1Føδ*O

Try it online! Footer formats output as a grid.

2āDδQ pushes the identity matrix to the stack:
2ā range from 1 to the length of the second input
DδQ equality table with itself

1F iterate first input times:
øδ*O multiply current matrix with input matrix:
ø transpose current matrix
δ* 3d multiplication table with input and transposed matrix
O sum along last axis