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]

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.

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

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";
[download]

bw, bliako

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)
#
# ...
[download]

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.

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

Cheers,
Rob