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; } } }

Download the source: qr_householder_bidiagonalization.cc.bz2