# Bidiagonalization using Householder transformations

The previous post was a discussion on employing Householder transformations to perform a QR decomposition. This post will be short. I've had this code lying around for a while now and thought I would make it available. The process of bidiagonalization using Householder transformations amounts to nothing more than alternating by left and right transformations.

The cMatrix::householderBidiagonalization() method:

void cMatrix::householderBidiagonalization(cMatrix& Q, cMatrix& R, cMatrix& S) {
double mag, alpha;

cMatrix u(m, 1),  v(m, 1),
u_(n, 1), v_(n, 1);

cMatrix P(m, m),  I(m, m),
P_(n, n), I_(n, n);

Q = cMatrix(m, m);
R = *this;
S = cMatrix(n, n);

for (int i = 0; i < n; i++) {
u.zero(); v.zero();

mag = 0.0;
for (int j = i; j < m; j++) {
u.A[j] = R.A[j * n + i];
mag += u.A[j] * u.A[j];
}
mag = sqrt(mag);

alpha = u.A[i] < 0 ? mag : -mag;

mag = 0.0;
for (int j = i; j < m; j++) {
v.A[j] = j == i ? u.A[j] + alpha : u.A[j];
mag += v.A[j] * v.A[j];
}
mag = sqrt(mag);

if (mag > 0.0000000001) {
for (int j = i; j < m; j++) v.A[j] /= mag;

P = I - (v * v.transpose()) * 2.0;

R = P * R;
Q = Q * P;
}
/////////////////////////
u_.zero(); v_.zero();

mag = 0.0;
for (int j = i + 1; j < n; j++) {
u_.A[j] = R.A[i * n + j];
mag += u_.A[j] * u_.A[j];
}

mag = sqrt(mag);

alpha = u_.A[i + 1] < 0 ? mag : -mag;

mag = 0.0;
for (int j = i + 1; j < n; j++) {
v_.A[j] = j == i + 1 ? u_.A[j] + alpha : u_.A[j];
mag += v_.A[j] * v_.A[j];
}
mag = sqrt(mag);

if (mag > 0.0000000001) {
for (int j = i + 1; j < n; j++) v_.A[j] /= mag;

P_ = I_ - (v_ * v_.transpose()) * 2.0;

R = R * P_;
S = P_ * S;
}
}
}