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

[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";
}

[bliako]

Hi, here is a modified version of your code to show the number of matrix multiplications too:

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

sub matrix_power {
    my ($M,$n,$nm) = @_;
    if($n<0) {
        return matrix_power(1/$M, -$n, $nm);
    } elsif ($n==0) {
        return 1;
    } elsif ($n==1) {
        return $M;
    } elsif (0 == $n % 2) {
    $$nm = $$nm + 1;
    print "(AxA)^".int($n/2)."\n";
        return matrix_power($M x $M, $n/2, $nm);
    } else {
    $$nm = $$nm + 2;
    print "Ax((AxA)^".int(($n-1)/2).")\n";
        return $M x matrix_power($M x $M, ($n-1)/2, $nm);
    }
}

my $x = pdl([0,1], [2,3]);
my $num_multiplications = 0;
my $power = 100;
print "M ** $power = " . matrix_power( $x, $power, \$num_multiplicatio
+ns )."\n";
print "NM=$num_multiplications\n";

bw, bliako

[vr]

Right, but couple of corrections/nitpicks, if I may:

if( $n < 0 ) {
    return matrix_power( inv( $M ), -$n );
} elsif ( $n == 0 ) {
    return identity( $M );
} elsif ( $n == 1 ) {
    return copy( $M );      # (1)
#
# ...

As to (1), I'm not sure about numpy, whether its function in question returns copy or clone in this case, or maybe they don't think it's important. But then, whole can of worms is there, e.g. data type of return piddle, what to do for singular matrices and negative exponent, and perhaps some others.

[pryrt]

Indeed. I did a rather naive interpretation of the wiki pseudocode. It would definitely need improvements (including those you showed) to be production-worthy.