gsl_linalg_complex_householder_transform

gsl_linalg_complex_householder_transform

[Mario Pernici]

[-- Attachment #1: Type: TEXT/PLAIN, Size: 1824 bytes --]


Dear Sirs,
 it seems that  gsl_linalg_complex_householder_transform
 does not follow Golub and Van Loan; is it done
 on purpose or there is a bug?

 Below I present the problem;
 in an attached files I give a modification of the
 code for gsl_linalg_complex_householder_transform
 according to Golub and Van Loan, the function
 gsl_linalg_complex_householder_hv
 and an example. Using the modified code for
 gsl_linalg_complex_householder_transform the example gives the
 expected results (es.1.out); the original code gives results
 not according to Golub and Van Loan (es.2.out).
 Best regards
   Mario

Problem:
 In the function
 double gsl_linalg_householder_transform(gsl_vector * v)
 there is a rescaling
  gsl_blas_dscal (1.0 / (alpha - beta), &x.vector);
 such that the first component of the householder vector is rescaled
 from (alpha - beta) to 1 (and then it is replaced by the norm).

 However this does not happen in
 gsl_complex
 gsl_linalg_complex_householder_transform (gsl_vector_complex * v).
 According to Golub and Van Loan, Problem P5.1.3, (their notation)
 householder matrix:
 P = 1 - beta u u^H
 P x = -e^{i theta} ||x||_2 * e_1
 householder vector:
 u = (x + e^{i theta} ||x||_2 * e_1 * s
 where x(1) = |x(1)| * e^{i theta}
 and where I added the scaling factor s, chosen to be
 s = 1/(x(1) + e^{i theta} ||x||_2) to set first component of the
 householder vector u(1) = 1.
 In the notation of linalg/householdercomplex.c the rescaling factor is
 s = 1/(alpha - e^{i theta} * beta_r);
 beta = (||x||_2 + |x(1)|)/||x||_2  becomes, in gsl notation,
 tau = (beta_r - |alpha|)/beta_r; notice that it is real (but I used
 a gsl_complex for it, not to change notation).

 In the gsl code the scaling factor is
 s = 1/(alpha - beta_r)
 so that u(1) is not rescaled to 1 when alpha = x(1) is not real.




[-- Attachment #2: Type: APPLICATION/x-gzip, Size: 1765 bytes --]

Re: gsl_linalg_complex_householder_transform

[Brian Gough]

Mario Pernici writes:
 > 
 > Dear Sirs,
 >  it seems that  gsl_linalg_complex_householder_transform
 >  does not follow Golub and Van Loan; is it done
 >  on purpose or there is a bug?

Hi,

I think I used the LAPACK method described in Lapack Working Note 72
"The Computation of Elementary Matrices", as it is supposed to have
better properties than the one in Golub and Van Loan.

I don't remember the details now -- I should have put some comments in
the code.

I will take another look at it to make sure it is ok.  Thanks for the
bug reports, keep them coming.

-- 
Brian

Re: gsl_linalg_complex_householder_transform

[Brian Gough]

Mario Pernici writes:
 >  Now my examples for the functions in the section Tridiagonal
 >  Decomposition of Hermitian Matrices work without modifications of
 >  the code of gsl.  The key point is that gsl_complex
 >  gsl_linalg_complex_householder_transform (gsl_vector_complex * v)
 >  returns a complex number tau such that (I - tau^* v v^H) x = -/+
 >  ||x|| e_1 (see Working Note 72) while I thought the relation was
 >  (I - tau v v^H) x = -/+ ||x|| e_1 Maybe this should be added in
 >  the comments to the source code of this function.

Yes, I should have commented it before... I've added some comments
now, and the function gsl_linalg_complex_householder_hv.

Brian

gsl_linalg_complex_householder_transform

[Mario Pernici]

Hi Brian,
< I think I used the LAPACK method described in Lapack Working Note 72

 thank you for the reference!
 Now my examples for the functions in the section
 Tridiagonal Decomposition of Hermitian Matrices
 work without modifications of the code of gsl.
 The key point is that
 gsl_complex
 gsl_linalg_complex_householder_transform (gsl_vector_complex * v)
 returns a complex number tau such that
 (I - tau^* v v^H) x = -/+ ||x|| e_1
 (see Working Note 72)
 while I thought the relation was (I - tau v v^H) x = -/+ ||x|| e_1
 Maybe this should be added in the comments to the source code of
 this function.

 And you could add to householdercomplex.c the function
gsl_linalg_complex_householder_hv in analogy with householder.c .
======================
int
gsl_linalg_complex_householder_hv(gsl_complex tau, const 
gsl_vector_complex * v, gsl_vector_complex *  w)
{
  size_t i;
  gsl_complex z;

  if (GSL_REAL(tau) == 0.0 && GSL_IMAG(tau) == 0.0)
      return GSL_SUCCESS;

  /* z = v^H w */
  gsl_blas_zdotc(v, w, &z);

  /* w = w - tau * (v^H w) * v   */
  z = gsl_complex_mul(z, tau);
  for (i = 0; i < v->size; i++)
  {
     gsl_complex vi, wi;
         vi = gsl_vector_complex_get(v,i);
         wi = gsl_vector_complex_get(w,i);
         vi = gsl_complex_mul(vi, z);
         wi = gsl_complex_sub(wi, vi);
         gsl_vector_complex_set(w, i, wi);
   }

}

=============================
Bye
 Mario

gsl_linalg_complex_householder_transform

[Mario Pernici]

[-- Attachment #1: Type: TEXT/PLAIN, Size: 2868 bytes --]


Dear Sirs,
 it seems that the modification of 
gsl_linalg_complex_householder_transform
 works fine; here I give a similar correction in gsl_linalg_hermtd_decomp.

 As in the previous mail,
 householder matrix:
 P = 1 - beta u u^H
 P x = -e^{i theta} ||x||_2 * e_1

 so that in gsl_linalg_hermtd_decomp
 the phase e^{i theta} must be taken into account when setting
 gsl_vector_complex_set (&v.vector, 0, ei);
 see attached file.

-------------------
 In the manual there are misprints:
 Tridiagonal Decomposition of Hermitian Matrices
A hermitian matrix A can be factorized by similarity transformations
 into the form,
 A = U T U^T
(becomes A = U T U^H)
 where U is an unitary matrix and T is a real symmetric tridiagonal 
matrix.
 ( T is a hermitian tridiagonal matrix).
---------------------

 Below I give an example in ruby-gsl; comments are added, so it should
 be understandable without running it.

require 'assert'
require "GSL"; include GSL; include GSL::Math
require 'gsl_matrix'; require 'gsl_linalg'

# complex numbers 0 and 1
zero = Complex.new2(0,0)
one = Complex.new2(1,0)

# initialize a 4x4 hermitian matrix
a = Matrix_complex.new([2,0,  3,2, 4,1, 5,2],
                       [3,-2, 2,0, 8,1, 4,4],
                       [4,-1, 8,-1,3,0, 3,2],
                       [5,-2, 4,-4,3,-2, 1,0])

# check that it is hermitian
assert a == a.h
# get the output matrix a1 and the vector tau using hermtd_decomp
a1, tau = a.hermtd_decomp

# In the lower part of a1 are stored the elements of the tridiagonal
# matrix and the householder vectors
#
#       [t             ]
#  a1 = [t    t        ]
#       [v_0  t   t    ]
#       [v_0  v_1 t  t ]

# householder vectors
v_0 = Vector_complex.new2(zero, one, a1.get(2,0),a1.get(3,0))
v_1 = Vector_complex.new2(zero,zero,one,a1.get(3,1))

# householder matrices h_i = I - tau_i v_i v_i^H
# where tau_i = 2/(v_i^H v_i)
h_0 = v_0.house
h_1 = v_1.house
# Q
q = h_0 * h_1
# Q is unitary
assert q * q.h =~ Matrix_complex.identity(4)

# Get the tridiagonal matrix t1 out of a1
# first put to zero the elements in the lower part where are stored
# the householder vectors
a1[2,0] = zero; a1[3,0] = zero; a1[3,1] = zero

# get the full hermitian tridiagonal matrix from its lower part
t1 = a1.lower(true) + a1.lower(false).h

# t1 = Q^H * a * Q within precision 1.0e-14
assert t1 =~ [q.h * a * q, 1.0e-14]
p t1
#[ (2, 0)  (-6.3911, -4.26073)  (0, 0)  (0, 0)
#   (-6.3911, 4.26073)  (9.89831, 0)  (-5.12537, -5.45784)  (0, 0)
#   (0, 0)  (-5.12537, 5.45784)  (-0.50685, 0)  (1.10104, -2.20513)
#    (0, 0)  (0, 0)  (1.10104, 2.20513)  (-3.39146, 0)  ]

# check that tau is related to the norms of the householder vectors
n0 = v_0.dznrm2
n1 = v_1.dznrm2
assert float_equal(2/(n0**2),tau[0].real, 1.0e-15)
assert float_equal(2/(n1**2),tau[1].real, 1.0e-15)
p tau
#[(1.4694, 0)  (0.632016, 0)  (-1.54304, 0.69892)  ]




[-- Attachment #2: Type: TEXT/PLAIN, Size: 2654 bytes --]

/* in linalg/hermtd.c */

int 
gsl_linalg_hermtd_decomp (gsl_matrix_complex * A, gsl_vector_complex * tau)  
{
  if (A->size1 != A->size2)
    {
      GSL_ERROR ("hermitian tridiagonal decomposition requires square matrix",
                 GSL_ENOTSQR);
    }
  else if (tau->size + 1 != A->size1)
    {
      GSL_ERROR ("size of tau must be (matrix size - 1)", GSL_EBADLEN);
    }
  else
    {
      const size_t N = A->size1;
      size_t i;
  
      const gsl_complex zero = gsl_complex_rect (0.0, 0.0);
      const gsl_complex one = gsl_complex_rect (1.0, 0.0);
      const gsl_complex neg_one = gsl_complex_rect (-1.0, 0.0);

      for (i = 0 ; i < N - 2; i++)
        {
          gsl_vector_complex_view c = gsl_matrix_complex_column (A, i);
          gsl_vector_complex_view v = gsl_vector_complex_subvector (&c.vector, i + 1, N - (i + 1));
          gsl_complex c0 = gsl_vector_complex_get(&c.vector, i+1);
          gsl_complex tau_i = gsl_linalg_complex_householder_transform (&v.vector);
	  /* get the phase of the first element */
	  c0 = gsl_complex_div_real(c0, gsl_complex_abs(c0));
          /* Apply the transformation H^T A H to the remaining columns */

          if ( !(GSL_REAL(tau_i) == 0.0 && GSL_IMAG(tau_i) == 0.0)) 
            {
              gsl_matrix_complex_view m = 
                gsl_matrix_complex_submatrix (A, i + 1, i + 1, 
                                              N - (i+1), N - (i+1));
              gsl_complex ei = gsl_vector_complex_get(&v.vector, 0);
              gsl_vector_complex_view x = gsl_vector_complex_subvector (tau, i, N-(i+1));
              gsl_vector_complex_set (&v.vector, 0, one);
              
              /* x = tau * A * v */
              gsl_blas_zhemv (CblasLower, tau_i, &m.matrix, &v.vector, zero, &x.vector);

              /* w = x - (1/2) tau * (x' * v) * v  */
              {
                gsl_complex xv, txv, alpha;
                gsl_blas_zdotc(&x.vector, &v.vector, &xv);
                txv = gsl_complex_mul(tau_i, xv);
                alpha = gsl_complex_mul_real(txv, -0.5);
                gsl_blas_zaxpy(alpha, &v.vector, &x.vector);
              }
              
              /* apply the transformation A = A - v w' - w v' */
              gsl_blas_zher2(CblasLower, neg_one, &v.vector, &x.vector, &m.matrix);

              /* multiply ei by the phase */
	      ei = gsl_complex_mul(ei, c0);
              gsl_vector_complex_set (&v.vector, 0, ei);
            }
          
          gsl_vector_complex_set (tau, i, tau_i);
        }
      
      return GSL_SUCCESS;
    }
}