From: David D. <dav...@gm...> - 2009-07-15 18:35:12
|
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 |
From: Stefan A. <at...@cs...> - 2009-07-16 13:23:46
|
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 > |