Advanced Search

Question: MatrixPower doesn't give the right answer

Question:MatrixPower doesn't give the right answer

mmcdara 6730 Maple
linearalgebra matrixpower
February 28 2023
2

I incidently(*) discovered that using MatrixPower to find the square root matrix of a relatively small matrix, returns wrong results, even for symetric positive definite matrices.

Trying to understand how MatrixPower proceeds using showstat is not that easy, but something I'm almost sure is that it doesn't ise an Eigenvectors decomposition.

Here is an example where I'm seeking for a matrix M such that C = M^2 and C is symetric definite positive of dimension 40.
(I didn't check if problems already occurs for smaller sizes)

(*)  I was very bored by the prohibitive execution time of MatrixPower and looked for an alternative.

The matrix M that MatrixPower generates is far from verifying C = M^2 (using the option 'shape'='symmetric' doesn't improve anything).
The naive alternative based on an Eigenvectors decomposition works perfectly well while being 200 times faster.

> 

restart:
> 

with(LinearAlgebra):
with(Statistics):
> 

f := proc(p)
  Sample(Uniform(0, 1), [100, p]):
  CorrelationMatrix(%):
end proc:
> 

C := f(40):
> 

M_MP := CodeTools:-Usage( MatrixPower(C, 1/2) ):

memory used=382.38MiB, alloc change=12.01MiB, cpu time=6.05s, real time=6.26s, gc time=386.00ms

  
> 

LE := proc(R)
  local L, E;
  L, E := LinearAlgebra:-Eigenvectors(evalf(R)):
  L    := map(Re, L):
  E    := map(Re, E):
  E . DiagonalMatrix(sqrt~(L)) . E^(-1)
end proc:
> 

IsDefinite(C, query = 'positive_definite');

M_LE := CodeTools:-Usage( LE(C) ):

true

  

memory used=0.90MiB, alloc change=0 bytes, cpu time=29.00ms, real time=10.00ms, gc time=0ns

  
> 

map(fnormal, M_MP . M_MP - C):
%[1..5, 1..5];

Matrix([[-493.2020830, -3096.801269, 3350.121141, -2972.190814, -6160.797678], [-2436.661317, -6415.706483, 208.1172241, -3001.930739, 1332.904551], [202.2560565, 413.4962022, 799.7833045, -179.6552377, -1817.784899], [-2994.580353, -9862.899715, 2543.110865, -5679.831865, -2552.349252], [-98.63574955, 2078.317993, -3118.445704, 2420.432714, 5905.315145]])

 (1)
> 

map(fnormal, M_LE . M_LE - C):
%[1..5, 1..5];

Matrix([[0., -0., -0., 0., 0.], [-0., 0., 0., 0., 0.], [-0., 0., -0., 0., -0.], [0., 0., -0., 0., -0.], [-0., 0., -0., -0., -0.]])

 (2)
> 

 

Download MatrixPower_bug.mw
Please Wait...