[Bug-gsl] linalg/bidiag.c

[Mario Pernici]

gsl version: 1.3
Red Hat 7.3
gcc version 2.96

Dear Sirs,
 There are errors in the comments in linalg/bidiag.c
See attached the proposed correction and an example.
I have also a remark.
 As the example shows, the orthogonal matrix U obtained from
gsl_linalg_bidiag_unpack is a M x N matrix (linalg/bidiag.c line 169)
satisfying U^T U = I, but not U U^T = I,
which is different from the M x M orthogonal householder matrix U
obtained from sl_linalg_bidiag_decomp;
V is the same in both cases; in both cases one gets a bidiagonal
matrix, but they differ (B has an extra row of zeroes).
 In both cases B = U^T A V and A = U B V^T .
 The gsl_linalg_bidiag_unpack results seem to be fine,
 but maybe it could be mentioned in the manual that
 the two U's and B's are different in the two functions.

Best Regards
 Mario Pernici

===================================
Proposed correction:
/* Factorise a matrix A into
 *
 * A = U B V^T
 *
 *........................................

 * The full matrix for U can be obtained as the product
 *
 *       U = U_0 U_1 .. U_(N-1)
 *
 * where
 *
 *       U_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [0,..,0,1, A(i+1,i), A(i+2,i), ... , A(M-1,i)]
 *
 * The full matrix for V can be obtained as the product
 *
 *       V = V_0 .. V_(N-3)
 *
 * where
 *
 *       V_i = (I - tau_i * v_i * v_i')
 *
 * and where v_i is a Householder vector
 *
 *       v_i = [0,..0,1, A(i,i+2), A(i,i+3), ... , A(i,N-1)]
 *
 * See Golub & Van Loan, "Matrix Computations" (3rd ed), Algorithm 5.4.2 
*/


==================================
Here is the example; I hope you don't mind it is in ruby-gsl!
I added comments so you do not need to run it.
---------------------------
require 'assert'
require "GSL"; include GSL; include GSL::Math
require 'gsl_matrix'; require 'gsl_linalg'
eps = 1.0e-15

# M x N matrix A with M = 4, N = 3; see Golub p252
A = Matrix.new([ 1, 2, 3],
               [ 4, 5, 6],
               [ 7, 8, 9],
               [10,11,12])
m = A.size1
n = A.size2
puts "bidiag_decomp"
# tm is the M x N matrix in which are stored the bidiagonal matrix B
# and the householder vectors
tm, tau_U, tau_V = A.bidiag_decomp
# householder vectors U_i, i = 0,.., N-1
u_0 = tm.col(0); u_0[0] = 1
u_1 = tm.col(1); u_1[0] = 0; u_1[1] = 1
u_2 = tm.col(2); u_2[0] = 0; u_2[1] = 0; u_2[2] = 1

# U = U_0 U_1 .. U_(N-1)
U = u_0.house * u_1.house * u_2.house

# householder vectors V_i, i = 0,.., N-3
v_0 = tm.row(0); v_0[0] = 0; v_0[1] = 1

# V = V_0 .. V_(N-3)
V = v_0.house

# B = U^T A V
B = U.t * A * V

# check that B is bidiagonal within precision 1.0e-14
(2..n-1).each{|i| assert B.superdiagonal(i) =~ [Vector::zero(n-i), 
1.0e-14]}
(1..m-1).each{|i| assert B.subdiagonal(i) =~ [Vector::zero(m-i), 1.0e-14]}

# check that tau_U(i) = 2.0/(u_i^T u_i)
tau_U1 = [u_0,u_1,u_2].map{|e| 2.0/(e*e)}
assert tau_U =~ Vector.new(tau_U1)

# check that tau_V(i) = 2.0/(v_i^T v_i) for i=0,..,N-3
# there is an extra tau_V(N-2) = 0
tau_V1 = [v_0].map{|e| 2.0/(e*e)}
tau_V1 << 0

assert tau_V =~ [Vector.new(tau_V1), 1.0e-15]

# B =
#[ -1.288e+01 2.188e+01 -4.892e-16
#  -3.331e-16 2.246e+00 -6.133e-01
#  1.110e-16 -9.375e-17 -4.376e-16
# -2.220e-16 7.011e-16 1.174e-16 ]
# in Golub the sign of B(0,0) is wrong

# V =
#[ 1.000e+00 0.000e+00 0.000e+00
#  0.000e+00 -6.670e-01 -7.451e-01
#  0.000e+00 -7.451e-01 6.670e-01 ]

# U =
#[ -7.762e-02 -8.331e-01 5.421e-01 -7.829e-02
#  -3.105e-01 -4.512e-01 -6.644e-01 5.084e-01
#  -5.433e-01 -6.942e-02 -2.974e-01 -7.820e-01
#  -7.762e-01 3.124e-01 4.198e-01 3.519e-01 ]


# There are misprints in the entries of U in Golub
U_Golub = Matrix.new([-0.0776, -0.833, 0.392, -0.383],
                     [-0.3110, -0.451, -0.238, 0.802],
                     [-0.5430, -0.069, 0.701, -0.457],
                     [-0.7760, 0.312, 0.547, 0.37])
# U_Golub.t * A * V =
#[ -1.288e+01 2.187e+01 -1.575e-04
#  -1.110e-16 2.246e+00 -6.131e-01
#   9.817e+00 -1.689e+01 5.946e-02
#   3.326e+00 -5.413e+00 -6.407e-02 ]
# is not bidiagonal

puts "bidiag_unpack"
vu, vv, vd, vs = tm.bidiag_unpack(tau_U, tau_V)
bd = vu.t * A * vv

# check that A = U B V^T using the output of bidiag_unpack
assert Matrix::float_equal(vu * bd * vv.t, A, 1.0e-14)

# bd is a bidiagonal matrix; vd and vu are its diagonal and first
# superdiagonal within precision 1.0e-14
assert bd.superdiagonal(0) =~ [vd, 1.0e-14]
assert bd.superdiagonal(1) =~ [vs, 1.0e-14]
#bd =
#[ -1.288e+01 2.188e+01 -4.892e-16
#  -7.772e-16 2.246e+00 -6.133e-01
#  3.331e-16 6.074e-16 -3.202e-16 ]

# vv is equal to V within precision 1.0e-15
assert V =~ vv

# vu is different from U
# vu is a 4 x 3 matrix; the first 3 columns of U are equal to those on vu
# however U is a 4 x 4 matrix
3.times {|i| assert vu.col(i) =~ U.col(i)}

# vu^T vu = I within precision 1.0e-15
assert vu.t * vu =~ Matrix.identity(3)

# vu * vu.t =
#[ 9.939e-01 3.981e-02 -6.122e-02 2.755e-02
#  3.981e-02 7.415e-01 3.976e-01 -1.789e-01
#  -6.122e-02 3.976e-01 3.885e-01 2.752e-01
#  2.755e-02 -1.789e-01 2.752e-01 8.762e-01 ]

# The first three rows of B and bd are equal within precision 1.0e-15
3.times {|i| assert B.row(i) =~ bd.row(i)}