Re: In PDL, how to raise a matrix $m to a power $v

[vr]

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 :)

[talexb]

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

does give

2  3
6 11

and multiplying it again produces

 6 11
22 39

(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.

[pryrt]

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

[syphilis]

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]
]


Cheers,
Rob

[vr]

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.

[pryrt]

Okay, then, implementing a full matrix_power() with PDL using "same algorithm as python numpy" (https://docs.scipy.org/doc/numpy-1.12.0/reference/generated/numpy.linalg.matrix_power.html), which says "the power is computed by repeated matrix squarings and matrix multiplications" -- which is then the matrix version of https://en.wikipedia.org/wiki/Exponentiation_by_squaring.

#!/usr/bin/perl
# [id://1229234]
use warnings;
use strict;
use PDL;

sub matrix_power {
    my ($M,$n) = @_;
    if($n<0) {
        return matrix_power(1/$M, -$n);
    } elsif ($n==0) {
        return 1;
    } elsif ($n==1) {
        return $M;
    } elsif (0 == $n % 2) {
        return matrix_power($M x $M, $n/2);
    } else {
        return $M x matrix_power($M x $M, ($n-1)/2);
    }
}

my $x = pdl([0,1], [2,3]);
for(0 .. 3) {
    print "M ** $_ = " . matrix_power( $x, $_ ) . "\n";
}

[syphilis]

Aaah ... I had checked whether perchance PDL overloaded '.' for matrix multiplication, but forgot about the possibility of overloading 'x'.
Odd choices, methinks.

Cheers,
Rob