A simple search on MetaCPAN finds Math::Matrix.

Looks like there is no such function. To clarify, we want, e.g., matrix

0 1
2 3

when squared, to become

2  3 
6 11

Searching CPAN, (link that talexb provided), Math::MatrixReal can do it, but inefficiently -- O(n) rather than O(log n). And it's definitely not PDL.

Indeed, as you say, very simple and reasonably efficient function can be written, utilizing exponentiation by squaring, maybe peeking at how numpy does it, but, of course, it's so simple, that the latter may not be necessary :)

This is really intriguing topic, and fondly brings back mathematical techniques that I probably haven't used in decades. What fun!

Squaring

0 1
2 3
[download]

does give

2  3
6 11
[download]

and multiplying it again produces

 6 11
22 39
[download]

(I'm sure there's a good reason why the top row of the cube has the same values as the bottom row of the square.) The module I suggested doesn't provide a power function, so I suppose the developer would have to jerry-rig something up that multiplies a matrix by itself n-1 times.

Alex / talexb / Toronto

Thanks PJ. We owe you so much. Groklaw -- RIP -- 2003 to 2013.

can't you just multiply the PDL matrix by itself?

can't you just multiply the PDL matrix by itself?

You can ... but it gives the result the OP doesn't want:

C:\>perl -MPDL -le "$x = pdl([0, 1], [2,3]);print $x; print $x * $x; p
+rint $x ** 2;"

[
 [0 1]
 [2 3]
]

[
 [0 1]
 [4 9]
]

[
 [0 1]
 [4 9]
]
[download]


Cheers,
Rob

I wrote "squared" only as an example. What if someone wants to raise to power of 100?

syphilis, for matrix multiplication, as opposed to this, PDL overloads the x operator.

The obvious shortcut here is to decompose the power p to a set of even-numbers whose sum equals p. For example for p=10, decompose 10 to 4,6 (just one of many). Then we do the same for each number (4 and 6).

Said decomposition (p=a+b+c+...) can be used to raise matrix, m, to power p as: m^p = m^a x m^b x m^c ... The benefit is that some calculations can be cached and re-used, thus minimising the number of matrix multiplications.

For example, m^14 = m^12 x m^2 First calculate m^2 = m x m and cache it (that's 1 multiplication). Then calculate m^4 = m^2 x m^2 and cache it (that's another one). Then calculate m^6 = m^4 x m^2 and cache it (thats' another one). Then calculate m^8 = m^4 x m^4 and cache it (thats' another one). Then calculate m^12 = m^8 x m^4 and cache it (thats' another one). Then calculate m^14 = m^12 x m^2 and cache it (thats' another one). ... Then multiply the cached results with 1 more multiplications. Done with 7 multiplications! (as opposed to p-1=13 required by the brute-force approach) Now, there is scope for even more optimisation using the decomposition m^14 = m^4 x m^4 x m^4 x m^2 which calculates m^14 using only 5 multiplications! So, there is a naive randomised algorithm which picks plans at random and choses the one with the fewer number of multiplications (without doing the actual multiplications) That's the "plan" which is then returned as a code-ref which is executed to yield the matrix to the asked power using PDL's matrix mult (x)

Edit: number of multiplications for some large power:

• p=50000 => ~20
• p=10000 => ~17
• p=1000 => ~14

Which is no O(log p) (Edit 2: it looks like my method is O(5log(p)) for these trials) as vr suggested (and I am not disputing). So there is scope for more optimisation. Come to think about it, I am not sure why I chose to decompose p to even-numbers!

Update 3: I haven't, at the time of first wrote this, searched any math sources for an algorithm to solve m^p but 24h later I did and sure there is Exponentiation by squaring https://en.wikipedia.org/wiki/Exponentiation_by_squaring and in simpler form https://simple.wikipedia.org/wiki/Exponentiation_by_squaring. So keep dividing the power by 2 (or simply follow the binary-equivalent of the power) and if the result is even do one thing, odd do another, recursively, just like the solution posted below Re^4: In PDL, how to raise a matrix $m to a power $v which yields exactly the same number of multiplications as mine (for 50000, 10000, 1000) and keeps its cache in the stack :) So, my method can be used if one does not want to use a recursive function. But it has a random element in it which I am sure, now in hindsight, eliminate and follow the odd/even/divide-by-2 method to find the best decomposition. Nice little problem.

bw, bliako

See mpow. On the plus side, it allows strange but probably amazing additional things -- fractional powers, or (reading source) power being not scalar but another matrix. On the minus side, see this comment. For familiar integer powers, it does O(n) loop (line 544). And it's buggy: if power is 0, same matrix is returned (as if power was 1). To fix, replace line 539 with $ret = identity($m); and start loop (line 544) from 0. And use PDL::MatrixOps qw/identity/; somewhere. But formally, yes, such function does exist.

Good catch! I've fixed this in https://github.com/PDLPorters/pdl-linearalgebra/commit/a2d680992c1cd4a46a3009b8ae0b2f3e6ebf601e. In fact, to aid broadcasting, I used (using two PDL::LinearAlgebra utility functions recently added):

my $di = $_[0]->dims_internal;
my @dims = $m->dims;
# ...
$ret = identity($dims[$di]);
$ret = $ret->r2C if $m->_is_complex;
[download]

This has now been released as 0.30.

Please do report problems like this at https://github.com/PDLPorters/pdl-linearalgebra/issues, we really do care! As an example, #10 (the wrapper functions didn't work on native complex, though the underlying LAPACK-bindings did) turned into a bit of a death-march because it revealed lots of copy-paste (~6800 lines ended up as ~3400), and very minimal tests. But I felt this module was important enough to justify the effort.