Re: expm bug fix

[Ross A. Lippert]

De-balancing the result by a matrix multiply and inversion
is inelegant.  If you finish your patch to CMatrix.cc then I'll
take what you did to expm and use dgebak to do the de-balancing
and get new patches out in a couple of days.

However, there is still the problem in dMatrix::solve() in which
if will return a [] when non-[] data is given to it.  That needs 
to be dealt with too, and it is a pretty fundamental problem.

-r

"A. Scottedward Hodel" wrote:
> 
> I wasn't on the octave-maintainers list until recently;
> I've sent the fix below to John Eaton but thought I'd
> distribute it to the list as well.
> 
> >Actually the bug in expm can be recreated easily by doing this:
> >
> >octave:15> y = triu(randn(2),1); y(2,1) = 1e-32;
> >octave:16> expm(y)
> >ans = [](0x0)
> >
> >This problem also occurs in version 2.0.14.
> >
> >
> 
> I've solved the expm problem. I've fixed the real expm
> case and will port and check the changes in the complex
> expm case next.
> 
> The problem lies in the use of AEPBALANCE to construct
> the balancing matrix.
> 
> AEPBALANCE calls fortran routine dgebal which returns
> balancing parameters (int) ilo, (int) ihi, and (double*) pscale.
> 
> dMatrix.cc: expm was modified to call dgebal and dgebak
>  directly instead of calling AEPBALANCE.  This gives
>  expm access to the balancing transformation through
>  variables ilo, ihi, and pscale.
> 
>  Routine dgebal computes m := inv(d)*m*d where d is
>         represented in variables ilo, ihi, and pscale.
> 
>  Routine dgebak can be used to directly compute the
>  matrix products:
>  d*a      (side = 'R')
>  a*inv(d) (side = 'L')
> 
>         Unfortunately, inverse scaling requires computation of
>         the matrix product
> 
>         retval = inv(d)*retval*d
> 
>         and so I chose to compute the inverse balancing via
>         matrix product.  d and inv(d) are computed from dgebak as
>         described above.
> 
> This is sloppy, since it's certainly possible to modify
> pscale, ilo, and ihi in order to be able to get
> the correct operations out of dgebak (which would be more
> efficient than a matrix multiply).  However, my time's
> pressed, so this is the solution I'm going with for now.
> 
> Listed below are
> (1) The (simplistic) test m-file I used to check my expm fix, and
> (2) the patch to dMatrix.cc and CMatrix.cc
> 
> Question: AEPBAL calls (*current_liboctave_error_handler)
> after the call to dgebal.  Is this necessary?  I did it
> in dMatrix.cc but not in CMatrix.cc.
> 
> (1) expmtest.m [taken from the user's email]
> y = triu(randn(2),1)
> y(2,1) = 1e-32
> z = expm(y)
> 
> y = y+1e-12*sqrt(-1)
> z = expm(y)
> 
> (2) Patch file to liboctave/dMatrix.cc and liboctave CMatrix.cc
> ===================================================================
> RCS file: CDiagMatrix.cc,v
> retrieving revision 1.1
> diff -u -r1.1 CDiagMatrix.cc
> ===================================================================
> RCS file: CMatrix.cc,v
> retrieving revision 1.1
> diff -u -r1.1 CMatrix.cc
> --- 1.1 1999/06/03 14:12:42
> +++ CMatrix.cc 1999/06/04 16:48:05
> @@ -21,6 +21,8 @@
> 
>  */
> 
> +#undef DEBUG
> +
>  #if defined (__GNUG__)
>  #pragma implementation
>  #endif
> @@ -33,6 +35,10 @@
> 
>  #include <iostream.h>
> 
> +#ifdef DEBUG
> +#include <iomanip.h>
> +#endif
> +
>  // XXX FIXME XXX
>  #ifdef HAVE_SYS_TYPES_H
>  #include <sys/types.h>
> @@ -59,6 +65,15 @@
> 
>  extern "C"
>  {
> +  int F77_FCN (zgebal, ZGEBAL) (const char*, const int&, Complex*,
> +                                const int&, int&, int&, double*, int&,
> +                                long, long);
> +
> +  int F77_FCN (dgebak, DGEBAK) (const char*, const char*, const int&,
> +                                const int&, const int&, double*,
> +                                const int&, double*, const int&,
> +                                int&, long, long);
> +
>    int F77_FCN (zgemm, ZGEMM) (const char*, const char*, const int&,
>           const int&, const int&, const Complex&,
>           const Complex*, const int&,
> @@ -1592,10 +1607,59 @@
>      m.elem (i, i) -= trshift;
> 
>    // Preconditioning step 2: eigenvalue balancing.
> +  // code follows development in AEPBAL
> 
> -  ComplexAEPBALANCE mbal (m, "B");
> -  m = mbal.balanced_matrix ();
> -  ComplexMatrix d = mbal.balancing_matrix ();
> +  #ifdef DEBUG
> +  // save a copy of m for comparison later
> +  ComplexMatrix m_copy = m;
> +  #endif
> +
> +  Complex *mp = m.fortran_vec();
> +  int ilo, ihi, info;
> +  char job = 'B';   // FIXME: should pass job as a parameter in expm
> +  Array<double> scale(nc);
> +  double *pscale = scale.fortran_vec();
> +
> +  F77_XFCN (zgebal, ZGEBAL, (&job, nc, mp, nc, ilo, ihi,
> +                             pscale, info, 1L, 1L));
> +
> +  // construct balancing matrices dmat, dinv
> +  Matrix dmat = Matrix(nc, nc, 0.0);
> +  Matrix dinv = Matrix(nc, nc, 0.0);
> +  int ii;
> +  for( ii = 0 ; ii < nc ; ii++)
> +    dmat(ii,ii) = dinv(ii,ii) = 1.0;
> +
> +  // pscale contains real data, so just use dgebak
> +  // dgebak, dside=R => dmat := D*dmat
> +  char dside = 'R';
> +  F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
> +                             dmat.fortran_vec(), nc, info, 1L, 1L));
> +
> +  // dgebak, dside=L => dinv := dinv*D^{-1}
> +  dside = 'L';
> +  F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
> +                             dinv.fortran_vec(), nc, info, 1L, 1L));
> +
> +  #ifdef DEBUG
> +  ComplexMatrix mchk = dmat*m*dinv;
> +  cout << "check: balanced m, m_copy, d*m*dinv - m_copy =" << endl;
> +  for( ii = 0; ii < nc ; ii++)
> +    for( int jj = 0 ; jj < nc ; jj++)
> +      cout
> + << setw(3) << ii
> + << setw(3) << jj
> + << "("
> + << setw(12) << setprecision(4) << mchk(ii,jj).real() << ", "
> + << setw(12) << setprecision(4) << mchk(ii,jj).imag() << ") "
> + << "("
> + << setw(12) << setprecision(4) << m_copy(ii,jj).real() << ", "
> + << setw(12) << setprecision(4) << m_copy(ii,jj).imag() << ") "
> + << "("
> + << setw(12) << setprecision(4) << (mchk(ii,jj)-m_copy(ii,jj)).real() << ", "
> + << setw(12) << setprecision(4) << (mchk(ii,jj)-m_copy(ii,jj)).imag() << ") "
> + << endl;
> +  #endif
> 
>    // Preconditioning step 3: scaling.
> 
> @@ -1659,14 +1723,10 @@
>      }
> 
>    // Reverse preconditioning step 2: inverse balancing.
> -  // XXX FIXME XXX -- should probably do this with Lapack calls
> -  // instead of a complete matrix inversion.
> -
> -  retval = retval.transpose ();
> -  d = d.transpose ();
> -  retval = retval * d;
> -  retval = d.solve (retval);
> -  retval = retval.transpose ();
> +  // FIXME: need to figure out how to get dgebak to do
> +  // dmat*retval*dinv instead of dinv*retval*dmat
> +  // This works for now
> +  retval = dmat*retval*dinv;
> 
>    // Reverse preconditioning step 1: fix trace normalization.
> 
> ===================================================================
> RCS file: boolMatrix.cc,v
> retrieving revision 1.1
> diff -u -r1.1 boolMatrix.cc
> ===================================================================
> RCS file: chMatrix.cc,v
> retrieving revision 1.1
> diff -u -r1.1 chMatrix.cc
> ===================================================================
> RCS file: dDiagMatrix.cc,v
> retrieving revision 1.1
> diff -u -r1.1 dDiagMatrix.cc
> ===================================================================
> RCS file: dMatrix.cc,v
> retrieving revision 1.1
> diff -u -r1.1 dMatrix.cc
> --- 1.1 1999/06/03 14:12:42
> +++ dMatrix.cc 1999/06/04 16:34:21
> @@ -21,6 +21,8 @@
> 
>  */
> 
> +#undef DEBUG
> +
>  #if defined (__GNUG__)
>  #pragma implementation
>  #endif
> @@ -33,6 +35,10 @@
> 
>  #include <iostream.h>
> 
> +#ifdef DEBUG
> +#include <iomanip.h>
> +#endif
> +
>  #include "byte-swap.h"
>  #include "dMatrix.h"
>  #include "dbleAEPBAL.h"
> @@ -54,6 +60,15 @@
> 
>  extern "C"
>  {
> +  int F77_FCN (dgebal, DGEBAL) (const char*, const int&, double*,
> +                                const int&, int&, int&, double*,
> +                                int&, long, long);
> +
> +  int F77_FCN (dgebak, DGEBAK) (const char*, const char*, const int&,
> +                                const int&, const int&, double*,
> +                                const int&, double*, const int&,
> +                                int&, long, long);
> +
>    int F77_FCN (dgemm, DGEMM) (const char*, const char*, const int&,
>           const int&, const int&, const double&,
>           const double*, const int&,
> @@ -1326,86 +1341,168 @@
>   m.elem (i, i) -= trshift;
>      }
> 
> -  // Preconditioning step 2: balancing.
> -
> -  AEPBALANCE mbal (m, "B");
> -  m = mbal.balanced_matrix ();
> -  Matrix d = mbal.balancing_matrix ();
> -
> -  // Preconditioning step 3: scaling.
> -
> -  ColumnVector work(nc);
> -  double inf_norm;
> -
> -  F77_FCN (xdlange, XDLANGE) ("I", nc, nc, m.fortran_vec (), nc,
> -         work.fortran_vec (), inf_norm);
> -
> -  int sqpow = (int) (inf_norm > 0.0
> -       ? (1.0 + log (inf_norm) / log (2.0))
> -       : 0.0);
> -
> -  // Check whether we need to square at all.
> -
> -  if (sqpow < 0)
> -    sqpow = 0;
> -
> -  if (sqpow > 0)
> -    {
> -      double scale_factor = 1.0;
> -      for (int i = 0; i < sqpow; i++)
> - scale_factor *= 2.0;
> -
> -      m = m / scale_factor;
> -    }
> -
> -  // npp, dpp: pade' approx polynomial matrices.
> +  // Preconditioning step 2: balancing; code follows development
> +  // in AEPBAL
> 
> -  Matrix npp (nc, nc, 0.0);
> -  Matrix dpp = npp;
> +  #ifdef DEBUG
> +  Matrix m_copy = m;
> +  cout << "dMatrix.cc: debugging code m_copy made " << endl;
> +  #endif
> +
> +  double *p_m = m.fortran_vec();
> +
> +  Array<double> scale(nc);
> +  double *pscale = scale.fortran_vec();
> +
> +  int info, ilo, ihi;
> +  char job = 'B';  // both scale and permute
> +
> +  int ii;
> +  #ifdef DEBUG
> +  int jj;
> +  cout << "Original m" << endl;
> +  for (ii=0 ; ii < nc ; ii++)
> +  {
> +    for (jj=0 ; jj < nc ; jj++)
> +      cout << setw(12) << setprecision(4) << m(ii,jj) << " ";
> +    cout << endl;
> +  }
> +  #endif
> 
> -  // Now powers a^8 ... a^1.
> -
> -  int minus_one_j = -1;
> -  for (int j = 7; j >= 0; j--)
> -    {
> -      npp = m * npp + m * padec[j];
> -      dpp = m * dpp + m * (minus_one_j * padec[j]);
> -      minus_one_j *= -1;
> -    }
> -
> -  // Zero power.
> -
> -  dpp = -dpp;
> -  for (int j = 0; j < nc; j++)
> -    {
> -      npp.elem (j, j) += 1.0;
> -      dpp.elem (j, j) += 1.0;
> -    }
> -
> -  // Compute pade approximation = inverse (dpp) * npp.
> -
> -  retval = dpp.solve (npp);
> -
> -  // Reverse preconditioning step 3: repeated squaring.
> -
> -  while (sqpow)
> -    {
> -      retval = retval * retval;
> -      sqpow--;
> +  F77_XFCN (dgebal, DGEBAL, (&job, nc, p_m, nc, ilo, ihi, pscale, info, 1L, 
> 1L));
> +  if (f77_exception_encountered)
> +    (*current_liboctave_error_handler) ("unrecoverable error in dgebal");
> +  else
> +  {
> +    #ifdef DEBUG
> +    cout << "Balanced m" << endl;
> +    for (ii=0 ; ii < nc ; ii++)
> +    {
> +      for (jj=0 ; jj < nc ; jj++)
> +        cout << setw(12) << setprecision(4) << m(ii,jj) << " ";
> +      cout << endl;
> +    }
> +    cout << "Check: trying to construct m_copy from scale, and m" << endl;
> +    Matrix bal_m_copy = m;
> +    #endif
> +
> +    // construct balancing matrices D, Dinv; initialize dmat, dinv to 
> identity
> +    Matrix dmat = Matrix(nc,nc,0.0);
> +    Matrix dinv = Matrix(nc,nc,0.0);
> +    for(ii=0 ; ii < nc ; ii++)
> +      dmat(ii,ii) = dinv(ii,ii) = 1.0;
> +
> +    // dgebak, dside=R => dmat := D*dmat
> +    char dside = 'R';
> +    F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
> +                               dmat.fortran_vec(), nc, info, 1L, 1L));
> +
> +    // dgebak, dside=L => dinv := dinv*D^{-1}
> +    dside = 'L';
> +    F77_XFCN (dgebak, DGEBAK, (&job, &dside, nc, ilo, ihi, pscale, nc,
> +                               dinv.fortran_vec(), nc, info, 1L, 1L));
> +
> +    #ifdef DEBUG
> +    cout << "d,dinv,d*dinv" << endl;
> +    Matrix d_dinv = dmat*dinv;
> +    for (ii=0 ; ii < nc ; ii++)
> +    {
> +      for (jj = 0; jj < nc ; jj++)
> +        cout
> + << setw(3) << ii
> + << setw(3) << jj
> + << setw(12) << setprecision(4) << dmat(ii,jj) << " "
> + << setw(12) << setprecision(4) << dinv(ii,jj) << " "
> + << setw(12) << setprecision(4) << d_dinv(ii,jj) << " "
> + << endl;
> +    }
> +
> +    bal_m_copy = dmat*bal_m_copy*dinv;
> +
> +    // print out difference
> +    cout << "dMatrix:expm: dinv*m*d - orig_m = " << endl;
> +    for (ii=0 ; ii < nc ; ii++)
> +    {
> +      for (jj = 0; jj < nc ; jj++)
> + cout << setw(3) << ii
> + << setw(3) << jj
> + << setw(12) << setprecision(4) << bal_m_copy(ii,jj) << " "
> + << setw(12) << setprecision(4) << m_copy(ii,jj) << " "
> + << setw(12) << setprecision(4) << bal_m_copy(ii,jj) - m_copy(ii,jj)
> + << endl;
>      }
> -
> -  // Reverse preconditioning step 2: inverse balancing.
> +    #endif
> 
> -  retval = retval.transpose();
> -  d = d.transpose ();
> -  retval = retval * d;
> -  retval = d.solve (retval);
> -  retval = retval.transpose ();
> -
> -  // Reverse preconditioning step 1: fix trace normalization.
> -
> -  if (trshift > 0.0)
> -    retval = exp (trshift) * retval;
> +    // Preconditioning step 3: scaling.
> +
> +    ColumnVector work(nc);
> +    double inf_norm;
> +
> +    F77_FCN (xdlange, XDLANGE) ("I", nc, nc, m.fortran_vec (), nc,
> +           work.fortran_vec (), inf_norm);
> +
> +    int sqpow = (int) (inf_norm > 0.0
> +         ? (1.0 + log (inf_norm) / log (2.0))
> +         : 0.0);
> +
> +    // Check whether we need to square at all.
> +
> +    if (sqpow < 0)
> +      sqpow = 0;
> +
> +    if (sqpow > 0)
> +      {
> +        double scale_factor = 1.0;
> +        for (int i = 0; i < sqpow; i++)
> +   scale_factor *= 2.0;
> +
> +        m = m / scale_factor;
> +      }
> +
> +    // npp, dpp: pade' approx polynomial matrices.
> +
> +    Matrix npp (nc, nc, 0.0);
> +    Matrix dpp = npp;
> +
> +    // Now powers a^8 ... a^1.
> +
> +    int minus_one_j = -1;
> +    for (int j = 7; j >= 0; j--)
> +      {
> +        npp = m * npp + m * padec[j];
> +        dpp = m * dpp + m * (minus_one_j * padec[j]);
> +        minus_one_j *= -1;
> +      }
> +
> +    // Zero power.
> +
> +    dpp = -dpp;
> +    for (int j = 0; j < nc; j++)
> +      {
> +        npp.elem (j, j) += 1.0;
> +        dpp.elem (j, j) += 1.0;
> +      }
> +
> +    // Compute pade approximation = inverse (dpp) * npp.
> +
> +    retval = dpp.solve (npp);
> +
> +    // Reverse preconditioning step 3: repeated squaring.
> +
> +    while (sqpow)
> +      {
> +        retval = retval * retval;
> +        sqpow--;
> +      }
> +
> +    // Reverse preconditioning step 2: inverse balancing.
> +    retval = dmat*retval*dinv;
> +
> +    // Reverse preconditioning step 1: fix trace normalization.
> +
> +    if (trshift > 0.0)
> +      retval = exp (trshift) * retval;
> +  } // else of dgebal exception test
> 
>    return retval;
>  }
> 
> A S Hodel Assoc. Prof. Dept Elect Eng, Auburn Univ,AL  36849-5201
> On leave at NASA Marshall Space Flight Center (256) 544-1426
> Address until 15 Mar 2000:Mail Code TD-55, MSFC, Alabama, 35812
> http://www.eng.auburn.edu/~scotte