Re: [Vxl-users] Matrix powers

[Stefan A. <at...@cs...>]

I can't find one, but the best route for the general case is through
the matrix exponential. This is far from trivial, but perhaps netlib
has a method for it that can be wrapped (/netlib/ieeecss/cascade seems
to have a Pade approximation for the exponential).

Look at http://www.cs.cornell.edu/cv/researchpdf/19ways+.pdf   , most
of which will apply to the matrix power as well - repeated squaring
and other tricks should apply.

The code below would work for diagonalizable matrices. The inverse at
the end is unnecessary, V.transpose() should be the same.



On Wed, Jul 15, 2009 at 1:35 PM, David Doria<dav...@gm...> wrote:
> Is there a function to raise a vnl_matrix to a power?
>
> We could use the diagonalization method:
>
>        vnl_matrix<double> MatrixPower(const vnl_matrix<double> &M, const
> double MatPow)
>        {
>                vnl_symmetric_eigensystem<double> Eigs(M);
>
>                vnl_matrix<double> V = Eigs.V;
>                vnl_diag_matrix<double> D = Eigs.D;
>
>                vnl_matrix<double> DRaised(D.rows(), D.columns(), 0.0);
>                for(unsigned int i = 0; i < D.rows(); i++)
>                {
>                        DRaised(i,i) = pow(D(i,i), MatPow);
>                }
>
>                vnl_matrix<double> Raised = V * DRaised * vnl_matrix_inverse<double>(V);
>
>                return Raised;
>        }
>
> The types get more complicated (i.e. complex numbers) for non-integer
> powers, but this works fine for integer powers.
>
> Thanks,
>
> David
>
> ------------------------------------------------------------------------------
> Enter the BlackBerry Developer Challenge
> This is your chance to win up to $100,000 in prizes! For a limited time,
> vendors submitting new applications to BlackBerry App World(TM) will have
> the opportunity to enter the BlackBerry Developer Challenge. See full prize
> details at: http://p.sf.net/sfu/Challenge
> _______________________________________________
> Vxl-users mailing list
> Vxl...@li...
> https://lists.sourceforge.net/lists/listinfo/vxl-users
>