[Octave-bug-tracker] [bug #60738] logm returning incorrect result with some real non-symmetric matrices

[anonymous]

Follow-up Comment #8, bug #60738 (project octave):

I have attached the logm.m file. I will ask you to revert the changes just
applied and follow the attached file. I am confused why the patch did not
work, I updated my local copy just before starting working. The way the patch
marked the differences would have been very hard to follow as I effectively
made two paths with an if statement, the original code was essentially
unmodified and I put in a second path for Hermitian Matrices.

The following will describe how the current algorithm works:

It performs the following integral, (a definition of the log function)

ln(A)=int_I^A y^-1 dy

It uses a quadrature with limits between 0 and 1 so if we apply a linear
coordinate transformation to change the lower limit from I to 0 and the upper
limit from A to 1 using y=(A-I)x+I we get the following integral

ln(A)=int_0^1 (A-I)((A-I)x+I)^-1 dx
The following lines in the function do this

s -= eye (n);
if (m > 1)
  s = logm_pade_pf (s, m);
endif

The following code is equivelent but much slower

n=3;
A=randn(n);
A*=A';#quad not a fan of complex numbers so make A symmetric and pos def
logmA=zeros(n);
for i=1:n
  for j=1:n
    logmA(i,j)=quad(@(x)((A-eye(n))/((A-eye(n)).*x+eye(n)))(i,j),0,1);
  end
end
logmA
max(abs(expm(logmA)-A)(:))

I also believe the following comment in the file is incorrect.

##LOGM_PADE_PF   Evaluate Pade approximant to matrix log by partial
fractions.
##   Y = LOGM_PADE_PF(A,M) evaluates the [M/M] Pade approximation to
##   LOG(EYE(SIZE(A))+A) using a partial fraction expansion.

I have shown above we are obtaining logm by its integration definition. The
Padé approximant is in the form (a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...)
where a and b are found by solving a system of equations formed from
log(x)=(a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...)
d(log(x))dx=d((a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...))dx
d^2(log(x))dx^2=d^2((a0+a1*x+a2*x^2+...)/(1+b1*x+b2*x^2+...))dx^2
...
where x is the location to form the approximant about x=I in our case.

If the method is obtaining the log with integration what is the while loop
doing with sqrtm? It is rescaling the matrix as square rooting a matrix makes
it closer to the identity allowing the integral to be found using fewer
quadrature points. As an extra bonus the first term of the Taylor series of
log ln(A+I)=sum_(k=1)^inf A^k/k is A so log(A) approx equal A-I when A is
almost the identity. Therefore, there is a small chance the rescaling loop may
also find the matrix log.

Why is the algorithm having issues with the log of ones(4)? There is a zero
divided by zero issue in sqrtm to fix change the algorithm for sqrtm from

A=ones(4);
n = rows (A);
[u T]=schur(A,'complex');
for j = 1:n
  T(j,j) = sqrt (T(j,j));
  for i = j-1:-1:1
    T(i,j) /= (T(i,i) + T(j,j));
    k = 1:i-1;
    T(k,j) -= T(k,i) * T(i,j);
  endfor
endfor
sqrtAA=u*T*u'
norm(sqrtAA^2-A,"fro")

to

A=ones(4);
n = rows (A);
[u T]=schur(A,'complex');
for j = 1:n
  T(j,j) = sqrt (T(j,j));
  for i = j-1:-1:1
    if T(i,j)~=0
      T(i,j) /= (T(i,i) + T(j,j));
    end
    k = 1:i-1;
    T(k,j) -= T(k,i) * T(i,j);
  endfor
endfor
sqrtAA=u*T*u'
norm(sqrtAA^2-A,"fro")

or change line 79 of sqrtm.cc from

const element_type colji = colj[i] /= (coli[i] + colj[j]);

to

if (colj[j] != zero){
  const element_type colji = colj[i] /= (coli[i] + colj[j]);
}

Also, ones(4) is a singular matrix so the inversion required in the
integration method is still going to cause issues.

Why is the alternative algorithm better for Hermitian matrixes?

The Schur decomposition A=U*S*U' is the first step of the existing algorithm.
S is normally upper triangular but it is diagonal for Hermitian Matrices,
therefore, when A is Hermitian (symmetric for real matrixes) the Schur
decomposition is equivalent to the eigendecomposition. Matrix functions are
easily applied to diagonal matrices logm(A)=U*diag(log(diag(S)))*U'. This only
applies when A is Hermitian and therefore the Schur decomposition results in S
being diagonal.




(file #51578)
    _______________________________________________________

Additional Item Attachment:

File name: logm.m                         Size:5 KB
    <https://file.savannah.gnu.org/file/logm.m?file_id=51578>



    _______________________________________________________

Reply to this item at:

  <https://savannah.gnu.org/bugs/?60738>

_______________________________________________
  Message sent via Savannah
  https://savannah.gnu.org/