Find the Matrix Power

MATL, 20 bytes

XJZyXyi:"!J2$1!*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×1×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 1×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
  2$1!  %   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

Jelly, 17 16 15 bytes

Z×'³S€€
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.

Z×'³S€€  Helper link. Argument: B (matrix)

Z        Zip; transpose rows and columns of B.
   ³     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.

APL, 32 31 chars

{⍺=0:(⍴⍵)⍴1⍨1+≢⍵⋄+.×⍨⍣(⍺-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.