Matrix raised to a matrix: $M^N$, is this possible? with $M,N\in M_n(\Bbb K).$ I was wondering if there is such a valid operation as raising a matrix to the power of a matrix, e.g. vaguely, if $M$ is a matrix, is
$$
M^N
$$
valid, or is there at least something similar?  Would it be the components of the matrix raised to each component of the matrix it's raised to, resulting in again, another matrix?
Thanks,
 A: It is possible to define the exponential
$$\exp(M) = \sum_{n \ge 0}^{\infty} \frac{M^n}{n!}$$
of any matrix using power series. Similarly, it is possible to define the logarithm
$$\log(I + M) = \sum_{n \ge 1} \frac{(-1)^{n-1} M^n}{n}$$
when the above series converges (this is guaranteed for example if the largest singular value of $M$ is less than $1$). We can therefore define
$$M^N = \exp(N \log M)$$
by imitating the identity $a^b = e^{b \log a}$ for, say, positive reals, but this won't have good properties unless $N$ and $M$ commute, I think. It's better to consider the exponential and logarithm separately. 
As I have discussed elsewhere on math.SE, the fact that the ordinary exponential takes two inputs which are the same type is misleading. Most (but not all) "exponential-type" operations in mathematics take two inputs which are different types. 
A: 
Please See first some quick explicit computations below:

I think the situation for to $2\times 2$ is more comprehensible. In oder to have a glimpse for the general case as mentioned by @Qiaochu Yuan
We know that the complex plan $\Bbb C$,
$$\Bbb C= \{a+ib: a,b\in\Bbb R\} $$ is isomorphic (as a field) to the to fields
$$G_2(\Bbb R) =\left\{\left(\begin{smallmatrix} a &b\\ -b&a\end{smallmatrix}\right): a,b\in\Bbb R\right\}$$
where one  merely identify
$$ 1 \equiv  \left(\begin{matrix}  1& 0\\0&1 \end{matrix}\right)~~~\text{and}~~~~ i \equiv  \left(\begin{matrix}  0&-1\\1&0 \end{matrix}\right)$$

Now notion of exponential makes senses for complex numbers up to branch cut of the logarithm. Namely is not surprising to write $$ z^w
~~\text{for}~~~z = a+ib, w =x+iy\in\Bbb C\setminus\{0\}.$$
which from the identification above readily represent the expression,
$$ \color{blue}{z^w\equiv \left(\begin{smallmatrix} a &-b\\ b&a\end{smallmatrix}\right)^\left(\begin{smallmatrix} x &-y\\ y&x\end{smallmatrix}\right)~\text{for}~~~a,b, x, y\in\Bbb R\setminus\{0\}}$$This expression is less abstract(but more particular and restricted) than, using the Taylor expansion for exponential and logarithm.


Another advantage of this expression is that, as opposed to the general cases where the expression of the $\log(M+I)$
$$\log(I + M) = \sum_{n \ge 1} \frac{(-1)^{n-1} M^n}{n}$$
converges only when $\|M\|< 1$, for the matrices taken in the fields $G_2(\Bbb R)$ such restriction is not required.

The above identification can be useful also to quickly compute the $n^{th}$ power of matrices(see here).
Application 1: Using the above identification we have,
$$ \displaystyle \left(\begin{smallmatrix} 0 &-1\\ 1&0\end{smallmatrix}\right)^\left(\begin{smallmatrix} 0 &-1\\ 1&0\end{smallmatrix}\right)\equiv i^i = e^{-\frac{\pi}{2}}\equiv \left(\begin{smallmatrix} e^{-\frac{\pi}{2}} &0\\ 0&e^{-\frac{\pi}{2}}\end{smallmatrix}\right) $$
Application 2: See here: :For each $a \in \mathbb{R}$ evaluate $ \lim\limits_{n \to \infty}\left(\begin{smallmatrix}1&\frac{a}{n}\\\frac{-a}{n}&1\end{smallmatrix}\right)^n$
$$\begin{align}\displaystyle \lim_{n \to \infty}\left(\begin{matrix} 1&\dfrac{a}{n}\\\dfrac{-a}{n}&1\end{matrix}\right)^{\left(\begin{smallmatrix} n &0\\ 0&n \end{smallmatrix}\right)} &\equiv \displaystyle \lim_{n \to \infty}\color{red}{\left(1-\dfrac{ai}{n}\right)^n}  \\&= \left(\begin{matrix} \cos a&\sin a\\-\sin a&\cos a\end{matrix}\right).\end{align}$$
Application 3:
$$ \displaystyle\left(\begin{smallmatrix} 1 &1\\ -1&1\end{smallmatrix}\right)^\left(\begin{smallmatrix} 0 &-1\\ 1&0 \end{smallmatrix}\right)\equiv (1+i)^i =  \left(\sqrt{2}e^{i\frac{\pi}{4}}\right)^i = (\sqrt{2})^ie^{-\frac{\pi}{4}}\\\equiv \displaystyle\left(\begin{smallmatrix} e^{-\frac{\pi}{4}} \cos\log \sqrt{2} &e^{-\frac{\pi}{4}}\sin\log \sqrt{2}\\ -e^{-\frac{\pi}{4}}\sin\log \sqrt{2}&e^{-\frac{\pi}{4}}\cos\log \sqrt{2}&\end{smallmatrix}\right) $$
Application4: with help of Application 3: compute
$$ \displaystyle\left(\begin{smallmatrix} 1 &1\\ -1&1\end{smallmatrix}\right)^\left(\begin{smallmatrix} 2 &1\\ -1&2 \end{smallmatrix}\right) $$
This last one is left as an exercise to make sure the OP and future readers understood the rule of computation going on.
A: In addition to the other answers, it is possible to rise any square matrix to the power of any other square matrix (even of different order).
I use the following Mathematica code for this:
Clear["Global`*"]
$PrePrint =.; FormatValues@Matrix = {};
Hold[MakeBoxes[Matrix[a_], StandardForm] ^:= 
    Block[{Internal`$ConvertForms = {}}, 
     TemplateBox[{MakeBoxes[MatrixForm@a][[1, 1, 3]]}, "Matrix", 
      DisplayFunction -> (RowBox@{style@"(", "\[NoBreak]", #, 
           "\[NoBreak]", style@")"} &)]]] /. 
  style[string_] -> StyleBox[string, FontColor -> Red] // ReleaseHold
Unprotect[Dot];
Dot[x_?NumberQ, y_] := x y;
Dot[A_?SquareMatrixQ, B_?SquareMatrixQ] := 
  Matrix[KroneckerProduct[A, 
      IdentityMatrix[LCM[Length[A], Length[B]]/Length[A]]] . 
     KroneckerProduct[B, 
      IdentityMatrix[LCM[Length[A], Length[B]]/Length[B]]]] /; 
   Length[A] != Length[B];
Protect[Dot];
Unprotect[Power];
Power[0, 0] = 1;
Protect[Power];
Matrix /: Matrix[x_?MatrixQ] := 
  First[First[x]] /; x == First[First[x]] IdentityMatrix[Length[x]];
Matrix /: NonCommutativeMultiply[Matrix[x_?MatrixQ], y_] := 
  Dot[Matrix[x], y];
Matrix /: NonCommutativeMultiply[y_, Matrix[x_?MatrixQ]] := 
  Dot[y, Matrix[x]];
Matrix /: Dot[Matrix[x_], Matrix[y_]] := Matrix[x . y];
Matrix /: Matrix[x_] + Matrix[y_] := Matrix[x + y];
Matrix /: x_?NumericQ + Matrix[y_] := 
  Matrix[x IdentityMatrix[Length[y]] + y];
Matrix /: x_?NumericQ Matrix[y_] := Matrix[x y];
Matrix /: Matrix[x_]*Matrix[y_] := Matrix[x . y] /; x . y == y . x;
Matrix /: Power[Matrix[x_?MatrixQ], y_?NumericQ] := 
  Matrix[MatrixPower[x, y]];
Matrix /: Power[Matrix[x_?MatrixQ], Matrix[y_?MatrixQ]] := 
  Exp[Log[Matrix[x]] . Matrix[y]];
Matrix /: Log[Matrix[x_?MatrixQ], Matrix[y_?MatrixQ]] := 
 Log[Matrix[x]]^-1 . Log[Matrix[y]]
Matrix /: Dot[Matrix[A_?SquareMatrixQ], Matrix[B_?SquareMatrixQ]] := 
 Matrix[KroneckerProduct[A, 
    IdentityMatrix[LCM[Length[A], Length[B]]/Length[A]]] . 
   KroneckerProduct[B, 
    IdentityMatrix[LCM[Length[A], Length[B]]/Length[B]]]]

$Post2 = 
  FullSimplify[# /. 
     f_[args1___?NumericQ, Matrix[mat_], args2___?NumericQ] :> 
      Matrix[MatrixFunction[f[args1, #, args2] &, mat]]] &;
$Post = Nest[$Post2, #, 3] /. Dot -> NonCommutativeMultiply &;

Now one can use something like this:
Matrix[ {
    {0, 1},
    {1, 0}
   } ] ^ Matrix[ {
    {0, 0, 1},
    {1, 0, 0},
    {0, 1, 0}
   } ]

And get the result. Note though, the computation takes long time and the result may be huge in volume and include Root objects.
