# Effective matrix power like algorithm

## First example

Suppose you want to calculate the 6th power of some matrix $A$. The brute force attempt of doing this is considering $$(((AA)A)A)A)A$$ which requires a total of 5 matrix multiplications. However, this can be improved on by smart grouping of the matrices, for example $$((AA)A)((AA)A)$$ requires only 3 multiplications, since when the value of say $AA$ is known from the left parenthesis, it does not have to be recomputed on the right hand side. Using this type of caching can save a lot of computational power (logarithmic scaling?).

## Second example

A lot more computation is necessary if we're not multiplying matrices but differential forms, since each multiplication requires the calculation of a whole lot of permutations in addition to the multiplication of the components. For example, the Liouville form in classical mechanics is $$\omega^n\propto\underbrace{\omega\wedge\omega\wedge\cdots\wedge\omega}_{n~\text{factors}}$$ where each $\omega$ is a 2-form.

## The question

Since this already smells like a good application of Mathematica's "x := x =" trickery, I would like to solve this problem in general: How do I divide the operations up in a smart way such that I have to do as little computational steps as possible, and how do I implement this in Mathematica? Note that matrices were only an example application (otherwise the answer would be MatrixPower of course), I sumbled upon that problem a couple of times in the past in very different situations, and had to resort to possibly very ineffective implementations.

## Mathematical statement

Take an associative algebra $\mathcal A$. How can $a^n$ with $a\in\mathcal A,~n\in\mathbb N$ be calculated using as little computations as possible?

-
Could you provide another example? I find it hard to generalize something to only one instance. – Mr.Wizard Feb 28 '12 at 23:07
you might find this mathematica-journal.com/issue/v5i4/columns/wagner/… interesting, although it is not the same problem – acl Feb 28 '12 at 23:19
Even the "simple" case of Project Euler Problem #122 is not so simple. To be clear are you looking for the best possible result, or merely a good one? For example, is the "binary" method from that page good enough? – Mr.Wizard Feb 28 '12 at 23:21
@Mr.Wizard Oh, I didn't remember it was directly a problem on PE. Unfortunately, I don't have the solution so I can't access the solutions forum there. Are there any nice insights hidden in there? The binary method seems nice, but since this question is partially academic anyway, I'd like to get some optimal solution as well. – David Feb 28 '12 at 23:26
@Artes Nest would be precisely the brute force attempt I am trying to avoid: Nest[#.a&, a, 5]$=((((aa)a)a)a)a$, which is two many computations. (pun pun) ;-) – David Feb 28 '12 at 23:29

Note: The method I describe below does not find the optimal solution (i.e. the minimal number of multiplications), but it usually does better than the binary approach, and more importantly: it is very easy to extend it and optimize it for special cases. We could easily pre-compute those exponents for which it does worse than the simple binary method to say up to 1000 and special case those.

According to the Wikipedia article on the topic:

... the determination of a shortest addition chain seems quite difficult: no efficient optimal methods are currently known for arbitrary exponents, and the related problem of finding a shortest addition chain for a given set of exponents has been proven NP-complete.

Therefore trying to compute the optimial solution within the function (as opposed to precomputing it and making a lookup table) defeats the purpose of using this for optimization.

I did a version of this with C++ template metaprogramming here.

This is a direct translation:

Let f be an associative function, then:

Clear[pow]
(* Syntax: pow[exponent][base] *)
pow[n_][a_] /; Mod[n, 2] == 0 := With[{t = pow[n/2][a]}, f[t, t]]
pow[n_][a_] /; Mod[n, 3] == 0 := With[{t = pow[n/3][a]}, f[f[t, t], t]]
pow[n_][a_] := f[pow[n - 1][a], a]
pow[1][a_] := a

Testing:

f[x_, y_] := (Print["x"]; x y)

Then pow[10][2] gives

During evaluation of In[389]:= x

During evaluation of In[389]:= x

During evaluation of In[389]:= x

During evaluation of In[389]:= x

Out[389]= 1024

i.e. it was done in four multiplications. You can think a bit about which is the best subdivision (to 3 or to 2) for different numbers. It can be manually verified that for exponents under 100, in the case of 33 and 69 (and exponents reducible to these such as $67 = 2\times 33 +1$) it's better not to subdivide into 3 first. We can easily special-case pow for these by adding the following to the beginning of pow's definition:

pow[33][a_] := f[pow[32][a], a]
pow[69][a_] := f[pow[68][a], a]

Similarly, in the case of e.g. 82 it's better not to subdivide into 2 because $82 - 1 = 3^4$, but use

pow[82][a_] := f[pow[81][a], a]

Another example is that in the case of 85 and 95 it's better to subdivide into 5 first:

pow[85][a_] := pow[13][ pow[5][a] ]

But most of these are just tweaks that will not make a huge difference.

I believe that the value of my implementation lies in the ease of its extension for special cases like these for small exponents.

For the adventurous with a lot of free time, see this paper on Pippenger's algorithm.

-
What's the reason for the division by 2 and 3? They look like primes to me, would adding a line for 5 improve the code? – David Feb 29 '12 at 4:41
@David This is not the optimal solution, but it's close enough to it to be practical. Are you looking for a method that will compute the optimal solution for any $n$, or are you looking for one that is practical to use (and fast)? Including division by 3 will usually improve the number of f-applications needed, but not always. You can verify manually that for 33, 69, 129, 141, 177, etc. it won't, and if needed you can add these as special cases. If you add division by 5, it will usually worsen the situation, but for e.g. 85 and 95 it will improve it. – Szabolcs Feb 29 '12 at 9:37
@David I am not at all sure that there is a simple enough way or closed formula to figure out the optimal solution. – Szabolcs Feb 29 '12 at 9:43
This solution would be easily generalized for larger primes (using Fold), but that makes it slower according to my testing. Adding 5 results in another step necessary, adding up to the 19th prime makes calculations take many many steps (although still being much better than brute force). – David Feb 29 '12 at 21:31
@David No, don't use Fold: f[pow[4][a],a] (assuming binary is implemented) is better than f[f[f[f[a, a], a], a], a] with my current solution. Also, bigger primes will not give better results int he majority of cases. The first case where 5 gives a better result is 85. – Szabolcs Feb 29 '12 at 21:40

Here's a functional implementation of the binary method. The Fold operates only Log[2,n] times, performing either one or two "multiplications" each time.

power[f_][a_, n_Integer] := With[{b = Reverse@IntegerDigits[n, 2]},
Last@Fold[
With[{p = f[#1[[1]], #1[[1]]]}, {p,
If[#2 == 1,
If[#1[[2]] === Null, p, f[p, #1[[2]]]], #1[[2]]]}] &,
{a, If[First@b == 1, a, Null]}, Rest@b]]

Here's an example, using a pure function which logs the multiplications it is asked to perform as well as actually doing them. Of course, f could be any associative dyadic operation, and a could be any suitable expression, not just a scalar as shown.

power[(Print["(", #1, "\[Times]", #2, ")"]; Times[#1, #2]) &][a, 6]

This will output $a^6$, but also print three lines $(a \times a)$, $(a^2 \times a^2)$ and $(a^4 \times a^2)$.

I've used Null as a placeholder to mean "no result yet". One could certainly add an argument to specify the actual identity element of f instead. Special-casing the least significant bit of n avoids doing an extra unnecessary squaring of the $a^{2^k}$ at the very end.

-
Binary is so passè, try the Zeckendorf representation. – rcollyer Feb 29 '12 at 2:23

The answer depends on the structure of the algebra.

I therefore cannot give a universal answer. Instead, I will illustrate the structural dependence by means of examples. In response to the request for implementation, they also illustrate one possible way to exploit Mathematica's support for symbolic computations. (I welcome comments suggesting clearer or more efficient approaches.)

Take, for instance, the second example of the question. Here is an abbreviated implementation of the wedge product of alternating bilinear forms, enough to get us going:

SetAttributes[Wedge, Flat];
Wedge[element[\[Omega]___], element[\[Eta]___]]  :=
With[{\[Phi] = Flatten[{\[Omega], \[Eta]}]},
signature[\[Phi]]  (element @@  Union[\[Phi]])
];
Wedge /: Wedge[Times[x_, \[Omega]_element], y_] := Times[x, Wedge[\[Omega], y]];
Wedge /: Wedge[x_, Times[y_, \[Omega]_element]] := Times[y, Wedge[x, \[Omega]]];
Wedge[\[Omega]_, 0] := 0;
Wedge[0, \[Omega]_] := 0;

It requires calculation of the signature of a permutation of arbitrary sortable objects, extended to have the value $0$ whenever there's a repetition. Here's a reasonably efficient implementation:

signature[p_List] :=
With[{c =  First[PermutationCycles[Ordering[p]] /. Cycles -> List]},
(-1)^(Length[Flatten[c]] - Length[c])
];
signature[p_List] /; Length[Union[p]] != Length[p] := 0;

More generally, in an arbitrary algebra the product will be computed by means of a matrix of structure constants specifying the product of any two elements in a given basis. The foregoing is a programmatically simple way to provide those structure constants for $\Lambda(V)$ for any finite-dimensional vector space $V$ (note that the dimension of $\Lambda(V^n)$ is $2^n$.) For efficiency (in low dimensions, anyway) we could compute these structure constants once and for all and cache them.

The trick to computing powers in algebras with lots of zeros and other symmetries in their structure matrices is to simplify as you go, as here:

Clear[WedgePower];
WedgePower[\[Omega]_, n_] /; n >= 2 :=
Nest[Distribute[Wedge[\[Omega], #], Plus, Wedge] &, \[Omega], n - 1];
WedgePower[\[Omega]_, 1] := \[Omega];
WedgePower[\[Omega]_, 0] := 1;

Distribute is implicitly doing the work of assembling pairs of basis elements to be multiplied and then collecting comparable terms, such as $dx_1\wedge dx_2 = -dx_2\wedge dx_1$, and summing their coefficients. As this is built in to Mathematica, we can expect it to be reasonably efficient.

Take, for instance, an arbitrary 2-form $\omega$ in $\Lambda(V^{10})$. It will have up to $\binom{10}{2} = 45$ nonzero coefficients; its square will have up to $\binom{10}{4} = 210$ nonzero coefficients; its cube, $210$; its fourth power, $45$; and its fifth power must be a multiple of a single form. Because you're working in an algebra of dimension $2^{10}=1024$, in principle each successive product requires $(2^{10})^2$, or over a million, calculations (although most of them will be products of zeros). Here, the bilinearity hugely reduces that number: the calculations needed to compute generic squares, cubes, fourth, and fifth powers of a 2-form are $2025$, $11475$, $20925$, and $22950$, respectively.

Trying to save time by breaking up powers into smaller groups (via the "binary method" or otherwise) can backfire. Continuing this example, computing $\omega^4$ via successive squaring requires first $2025$ operations to compute $\omega^2$ and then $44100$ to square that, a total of $46125$ products. Compare this to the value of $20925$ needed for successive wedging by $\omega$ itself: it will take more than twice as much effort.

By choosing an appropriate basis, the Liouville form can be written in a special way as $p_1 \wedge q_1 + \cdots + p_k \wedge q_k$. This reduces the calculations still further. With $k=5$ (working once more in a ten-dimensional space), for instance, the powers require only 25, 50, 50, and 25 wedge operations, respectively. This suggests yet another line of attack: before computing powers of an element of an algebra, change the basis if possible to simplify the calculation. Once again, though, how this is done will depend on the algebraic structure.

-