5
$\begingroup$

I am working on a problem that involves taking fractional powers of particular matrices.

For the matrix A with 2 on the main diagonal and -1 on the sub and super diagonal (the finite difference matrix for the Dirichlet Laplacian in 1D), there was basically no difference between using MATLAB's sqrtm(A) function and diagonalizing the matrix ([V,D]=eig(A), A = V\D*V) and computing V\D^(1/2)*V even if the matrix was as large as 1000x1000. This made me think I could compute other fractional powers V\D^(1/n)*V that MATLAB does not have built in functions for.

However when I consider the matrix A for the Dirichlet Laplacian in 2D, there is a significant difference between MATLAB's sqrtm function and the diagonalization procedure I used above even if the matrix was as small as 100x100. In particular, MATLAB's sqrtm function is better and my results don't make any sense if I do it the other way.

I am wondering if anyone can explain why this happens and if there is a way to fix it so that I can compute fractional powers of a large matrix in MATLAB.

Improve this question
$\endgroup$
2
  • $\begingroup$ You can directly use something like A^(2/3). Is there a reason this is undesirable? $\endgroup$
    – Doug Lipinski
    Nov 19, 2014 at 12:15
  • $\begingroup$ It is not undesirable, I just didn't know it could be this easy! Someone has also informed me that expm(1/n*logm(A)) works too. Both of these are giving me matrices whose nth power differs from A by basically machine precision. Thanks again, if you submit that as an answer I will accept it. $\endgroup$
    – RHP
    Nov 19, 2014 at 17:42

1 Answer 1

Reset to default
4
$\begingroup$

MATLAB has the ability to compute fractional powers of matrices by using the normal ^ operator which calls the mpower function when either of the arguments is a matrix. You can read a discussion of when fractional powers of matrices can be computed here.

Your original issues may stem from the fact that you have mixed up the formula to compute matrices raised to a fractional power. The matrix should be diagonalized as $A = PDP^{-1}$. Then $A^\alpha = PD^\alpha P^{-1}$.

In MATLAB this can be done as follows (your code is a bit off):

[V,D]=eig(A);
Aalph = V*(D.^alph)/V;

Based on a quick test, it seems that this is exactly the method used by MATLAB inside mpower. The following code returns exactly 0 (no differences due to finite precision) on my machine.

A = rand(1000);
alph = 2/3;
[V,D]=eig(A);
Aalph = V*(D.^alph)/V;
B = A^(alph);
max(abs(Aalph(:)-B(:)))
Improve this answer
$\endgroup$
1
  • $\begingroup$ Thanks! This is what happens when you think you remember a formula and don't stop to check that it is correct. At least now I won't be forgetting V*D/V any time soon. Thanks again for clearing that up. $\endgroup$
    – RHP
    Nov 19, 2014 at 21:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.