Fractional Powers of a Matrix

[bibhukalyana]

Hi everyone ,
I want to calculate fractional power of a matrix. How i will do it ? Is there any class related to matrix in qt ?

with regards
bibhu

[Cruz]

Fractional power of a matrix? Is that element wise? Otherwise I can't imagine how that should work. Anyhow, Qt doesn't have anyhting matrix related, but Armadillo (http://arma.sourceforge.net/) is an excellent and simple matrix library.

[bibhukalyana]

Thanks for reply,

If A is a matrix of 2x2 then A^2 means multiply matrix A with itself (means A X A).

[Cruz]

Sure, but what would be A^1.5? Since you're saying fractional powers.

[bibhukalyana]

Actually i also do not know how to calculate A^1.5 .
It may be like A^(3/2) = sqrt(A^3).

Actually i want to apply a equation like g(x,y) = iFFT { FFT(u,v) X |FFT(u,v)|^k } where 0 < k <1 |FFT(u,v)| is also a u X v matrix.

[Cruz]

Oh in that case Armadillo is no help, it's only good for simple linear algebra. With fourier transformation I'm afraid you gonna have to try a heavier gun, such as the GSL (http://www.gnu.org/software/gsl/) library. It has matrix classes and readily implemented fourier transformation functions.

[stampede]

Actually i also do not know how to calculate A^1.5 .
It may be like A^(3/2) = sqrt(A^3).
Really ? :)


A^(3/2) = sqrt(A^3) // your equation
let B = A^3, then
sqrt(A^3) = sqrt(B) = B^(1/2) // and what could that be ? ;)

If I remember correctly (forgive me any mistakes, I've passed the linear algebra course 4 yeas ago), fractional power of a (square) matrix is computed this way:
Let A be n x n matrix, then eigenvalues of A are the roots of a polynomial P(x) = det(A-x*I), where I is identity matrix of the same size as A. For each eigenvalue there is a corresponding eigenvector v_i, such that:
(1) A * v_i = e_i * v_i (e_i - i-th eigenvalue)
Above equation applies to any Natural number n:
(2) A^n * v_i = (e_i ^ n) * v_i (this can be easily proven by induction)
We can generalize this to any Rational number r instead of n (I dont remember the proof).
So in order to calculate A^x, where x is Rational number you need to find just one eigenvalue and corresponding eigenvector.

Simple example:
A = ( 3 1, 2 2 ) (row1, row2) (ok, this example is prepared such that roots are easily computed, but its only to show the method)
P(x) = det( A - x*I ) = det( 3-x 1, 2 2-x ) = (3-x)(2-x) - 2 = x^2 - 5x + 4 = (x-4)(x-1)
P(x) == 0 <=> x in { 4, 1 }
Eigenvector v=(x,y) for value 4 (equation (1)):
A* (x, y) = 4*(x,y)
(3x + y, 2x + 2y) = (4x, 4y) <=> x=y, so we may take any value, lets use x=y=2
v = (2,2)
Lets compute A^(3/2) (equation (2)):
A^(3/2) * v = 4^(3/2) * v
A^(3/2) * (2,2) = (2*4^(3/2), 2*4^(3/2))
Let A^(3/2) = B = (a b, c d) =>
(a b, c d)*(2,2) = (2*4^(3/2),2*4^(3/2)) => (calculated on paper) =>
a = 5/3, b = 1/3, c = 1/3, d = 4/3 (if im not wrong)

So A^(3/2) = (5/3 1/3, 1/3 4/3)

I used to do the calculations much quicker, I'm getting old :P

[bibhukalyana]

Thanks stampede,
In my case the matrix A is a image of 256X256.So if i want to calculate the eigenvalue the i think i have to solve a polynomial of 256 power.How i will solve this ? According to my knowledge solving polynomial is slower.Can you tell me how to do this things ?

[stampede]

I don't think you should implement fractional matrix powers yourself if you care about speed. I'd just use one of available scientific libraries (like already linked GSL). My example was just to show the relevant math ( because you said "i also do not know how to calculate A^1.5", so now you know :) ).

[bibhukalyana]

I had used matlab to calculate the fractional power of a square matrix of 256 X 256.It is faster.So i thought there should be any algorithm to do this.Is there any alternative or any shortcut for this ?

[stampede]

It is faster
It is faster than what ? GSL ?

So i thought there should be any algorithm to do this
Yes, there is, in GSL (or alike) library.
There is no such algorithms like Fourier Transforms in Qt, remember that this is mostly a GUI framework, you need to mix Qt with external libraries to provide additional functionality (like GSL for advanced math, libraries for serial port communication, video processing etc... )

[SixDegrees]

Just a side note: GSL is offered under the Gnu GPL. Using it in your code poses a significant threat to retaining control of your code. There are other mathematical libraries out there that provide matrix capabilities. If you're planning on distributing your software, you'll probably want to consult an attorney before committing to the GSL, or to any other packaged software that may place restrictions on distribution.

[MarekR22]

Form mathematical point of view it is possible.
There is a way to apply any analytical function on matrix.
First you have to expand function to taylor-polynomial (http://en.wikipedia.org/wiki/Taylor%27s_theorem). Since you can't calculate fractional power of negative value central point must be some bigger value (1 is convinient)
From that point it is ease to use this polynomial to calculate value of function for matrix.
so assuming that power is >0 it will be something like that:

typedef SommatrixClas Matrix;

Matrix power(cons Matrix &A, qreal pow) {
Q_ASSERT(pow>0);
Q_ASSERT(A.columnCount()==A.rowCount());
Matrix result = A;
Matrix b = 1;
Matrix old = 0;
for (int i=1; !qFuzzyCompare(result, old); ++i) {
b *= A * (pow/i);
pow-= 1.0;
old = result;
result += pow;
}
return result;
}

[bibhukalyana]

thanks MarekR22

But what is the use of b.

[MarekR22]

Ha, your question coused that i fond error in my code :).

typedef SommatrixClas Matrix;

Matrix power(cons Matrix &A, qreal pow) {
Q_ASSERT(pow>0);
Q_ASSERT(A.columnCount()==A.rowCount());
Matrix result = A;
Matrix x = A - 1; // here A has to be "moved" to take into account x0!
Matrix b = 1;
Matrix old = 0;
for (int i=1; !qFuzzyCompare(result, old); ++i) {
b *= x * (pow/i);
pow-= 1.0;
old = result;
result += b; // here was mistake too
}
return result;
}

This is just a concept demonstration, for example qFuzzyCompare should be replaced with some better condition of ending polynomial calculation, implementation can also change depending on used matrix class.

[bibhukalyana]

Can you explain me what is the meaning of A - 1 and X * (pow/i) [ if A = {1 2, 3 4} then A - 1 = ? ] ?