Re: [ublas] matrix exponential

Subject: Re: [ublas] matrix exponential
From: Paul Leopardi (paul.leopardi_at_[hidden])
Date: 2010-10-26 18:38:18

On Thursday 21 October 2010 20:32:43 Kraus Philipp wrote:
> Am 21.10.2010 um 09:31 schrieb sguazt:
> > On Thu, Oct 21, 2010 at 8:52 AM, David Bellot
> > <david.bellot_at_[hidden]> wrote:
> > 
> > Actually, MATLAB uses Pad� approximation + scaling&squaring method
> > (http://www.mathworks.com/help/techdoc/ref/expm.html)
> > 
> > Also OCTAVE seems to use the same algorithm:
> > http://hg.savannah.gnu.org/hgweb/octave/file/304b0ed4ca56/scripts/linear-
> > algebra/expm.m
> 
> It seems the Pad approximation is the fastest one, so I would like to
> use them

See my GluCat code ( http://glucat.sf.net )
The latest version (0.5.1) uses the scaling and squaring Pade approximation
method with UBLAS matrices to calclate the exponential of a multivector. 
See function exp() in matrix_multi_imp.h (code reproduced below).
The Clifford algebra notation complicates thing somewhat, but you can get a 
good idea of the algorithm from reading the code. See also the references

[GL]: Gene H. Golub and Charles F. van Loan,
"Matrix Computations", 3rd ed., Johns Hopkins UP, 1996.
[H]: Nicholas J. Higham
"The Scaling and Squaring Method for the Matrix Exponential Revisited",
SIAM Journal on Matrix Analysis and Applications,
Vol. 26,  Issue 4 (2005), pp. 1179-1193.

Best, Paul

  /// Exponential of multivector
  template< typename Scalar_T, const index_t LO, const index_t HI >
  const matrix_multi<Scalar_T,LO,HI>
  exp(const matrix_multi<Scalar_T,LO,HI>& val)
  {
    // Scaling and squaring Pade' approximation of matrix exponential
    // Reference: [GL], Section 11.3, p572-576
    // Reference: [H]

    typedef numeric_traits<Scalar_T> traits_t;

    if (val.isnan())
      return traits_t::NaN();

    const Scalar_T scalar_val = scalar(val);
    const Scalar_T scalar_exp = traits_t::exp(scalar_val);
    if (traits_t::isNaN_or_isInf(scalar_exp))
      return traits_t::NaN();
    if (val == scalar_val)
      return scalar_exp;

    typedef matrix_multi<Scalar_T,LO,HI> multivector_t;
    multivector_t A = val - scalar_val;
    const Scalar_T pure_scale = A.norm();

    if (traits_t::isNaN_or_isInf(pure_scale))
      return traits_t::NaN();
    if (pure_scale == Scalar_T(0))
      return scalar_exp;

    const int ilog2_scale =
      std::max(0, traits_t::to_int(ceil((log2(pure_scale) +
                  Scalar_T(A.frame().count()))/Scalar_T(2)) - 3));
    const Scalar_T i_scale = traits_t::pow(Scalar_T(2), ilog2_scale);
    if (traits_t::isNaN_or_isInf(i_scale))
      return traits_t::NaN();

    A /= i_scale;
    multivector_t pure_exp;
    {
      const int q = 13;
      static Scalar_T c[q+1];
      if (c[0] != Scalar_T(1))
      {
        c[0] = Scalar_T(1);
        for (int
            k = 0;
            k != q;
            ++k)
          c[k+1] = c[k]*(q-k) / ((2*q-k)*(k+1));
      }
      const multivector_t& A2 = A*A;
      const multivector_t& A4 = A2*A2;
      const multivector_t& A6 = A4*A2;
      const multivector_t& U =     
            c[0]+A2*c[2]+A4*c[4]+A6*(c[6]+A2*c[8]+A4*c[10]+A6*c[12]);
      const multivector_t& AV =
         A*(c[1]+A2*c[3]+A4*c[5]+A6*(c[7]+A2*c[9]+A4*c[11]+A6*c[13]));
      pure_exp = (U+AV) / (U-AV);
    }
    for (int
        k = 0;
        k != ilog2_scale;
        ++k)
      pure_exp *= pure_exp;
    return pure_exp * scalar_exp;
  }