This article presents the implementation of QR decomposition in D. Some background implementation details will be left out, for example implementation of the matrix structure used, and elements of common linear algebra algorithms. The gemm BLAS algorithm will be called from the C Openblas library.
References consulted for this article are:
The article focuses on the implementation of the Gram-Schmidt, its modified from, and the Householder QR decomposition, the referred texts give a good theoretical and intuitive account for the QR decomposition of a matrix \(A = QR\).
The map functions I use here are not the same as those in the D standard library, but the intuition for their implementation is covered in the account for implementing the much more interesting mapply function in another article.
Otherwise I import the items:
import std.stdio: writeln;
import core.stdc.math: abs = fabsf, abs = fabs, abs = fabsl,
pow, pow = powf, pow = powl,
sqrt, sqrt = sqrtf, sqrt = sqrtl;
This is a function to calculate the norm of an array:
T norm(bool doSqrt = true, T)(in T[] x)
if(isFloatingPoint!(T))
{
T value = 0;
for(long i = 0; i < x.length; ++i)
{
value += x[i]*x[i];
}
static if(doSqrt)
return sqrt(value);
else
return value;
}
the in qualification on the argument ensures that the array x is not modified. Implementation of a dot product is given here:
T dot(T)(in T[] x, in T[] y)
if(isFloatingPoint!(T))
{
T value = 0;
for(long i = 0; i < x.length; ++i)
{
value += x[i]*y[i];
}
return value;
}
the function below returns whether all the elements in the array are zero, it uses an epsilon defined by \(T_{\epsilon}^{0.5}\) where \(T_{\epsilon}\) is the machine precision of the floating point type.
bool allZero(T)(T[] x)
{
enum T eps = T.epsilon^^cast(T)(0.5);
for(long i = 0; i < x.length; ++i)
{
if(abs(x[i]) < eps)
return true;
else
continue;
}
return false;
}
the function calls will be distinquished with:
enum QRAlgorithm {GramSchmidt, ModifiedGramSchmidt, Householder}
alias GramSchmidt = QRAlgorithm.GramSchmidt;
alias ModifiedGramSchmidt = QRAlgorithm.ModifiedGramSchmidt;
alias Householder = QRAlgorithm.Householder;
In the matrix \(A\) with number rows and columns given by \(m \times n\), the columns of the matrix can be represented by \(n\) vectors \(a_{1}, \ldots, a_{n}\). The Gram-Schmidt algorithms calculates \(q_{i}\), which are columns of the Q matrix where \(i = 1, \ldots, n\). The first value \(\tilde{q_{1}} = a_{1}\). The algorithm is given by:
for \(i = 1, \ldots, n\):
\(\tilde{q_{i}} = a_{i} - (q_{1}^{T}a_{i})q_{1} - \ldots - (q_{i - 1}^{T}a_{i})q_{i - 1}\)
If \(\tilde{q_{i}} = 0\) quit.
\(q_{i} = \tilde{q_{i}}/\lVert\tilde{q_{i}}\rVert\)
the \(R\) matrix is calculated while \(Q\) is being calculated. \(R_{ij} = q_{i}^{T}a_{j}\) and \(R_{ii} = \lVert\tilde{q_{i}}\rVert\)
The code for the QR algorithm is given by:
auto qr(QRAlgorithm algo, T)(Matrix!(T) A)
if(isFloatingPoint!(T) && (algo == GramSchmidt))
{
auto Q = Matrix!(T)(A.nrow, A.ncol);
auto R = zeros!(T)(A.ncol, A.ncol);
T _norm = norm(A[0]);
R[0, 0] = _norm;
Q[0] = map!((a) => a/_norm)(A[0]);
long n = A.ncol;
for(long i = 1; i < n; ++i)
{
T[] qi = A[i].dup;
T[] ai = A[i];
for(long j = 0; j < i; ++j)
{
auto rji = dot(Q[j], ai);
R[j, i] = rji;
auto tmp = map!((a) => a*rji)(Q[j]);
qi[] -= tmp[];
}
_norm = norm(qi);
R[i, i] = _norm;
Q[i] = map!((a) => a/_norm)(qi);
if(allZero(Q[i]))
{
assert(0, "Error, qi = 0");
}
}
return [Q, R];
}
we will run all three QR algorithms on the following matrix:
auto X = Matrix!(double)([0.8594598509, 0.8886035203, 0.8149294811, 0.7431045200, 0.8032585254,
0.0587533356, 0.7245921139, 0.5380305406, 0.7342256338, 0.6982547215,
0.7176400044, 0.0539911194, 0.3670289037, 0.9701228316, 0.8404100032,
0.4112932913, 0.3075223914, 0.5798244230, 0.0015286701, 0.7890766996,
0.9781337455, 0.2921431712, 0.0432923459, 0.9428416709, 0.9646959945,
0.0354323143, 0.4898468039, 0.4513681016, 0.2107982126, 0.4445287671,
0.8115565467, 0.7058405790, 0.5527189195, 0.5410537042, 0.9117912347,
0.1149175267, 0.8406228190, 0.6040554044, 0.4260203703, 0.2376075180,
0.2164094832, 0.1800869710, 0.7479251262, 0.0009715103, 0.8810979640,
0.8647838791, 0.5856765260, 0.0127644690, 0.5744975219, 0.1985024847], 10, 5);
writeln("X: ", X);
output:
0.8595 0.7176 0.9781 0.8116 0.2164
0.8886 0.05399 0.2921 0.7058 0.1801
0.8149 0.367 0.04329 0.5527 0.7479
0.7431 0.9701 0.9428 0.5411 0.0009715
0.8033 0.8404 0.9647 0.9118 0.8811
0.05875 0.4113 0.03543 0.1149 0.8648
0.7246 0.3075 0.4898 0.8406 0.5857
0.538 0.5798 0.4514 0.6041 0.01276
0.7342 0.001529 0.2108 0.426 0.5745
0.6983 0.7891 0.4445 0.2376 0.1985
The code to execute the Gram-Schmidt QR algorithm is given by:
writeln("qr!(GramSchmidt)(X): \n", qr!(GramSchmidt)(X));
the output obtained:
[Matrix(10 x 5)
0.3757 0.1337 0.4163 -0.07128 -0.02322
0.3884 -0.4844 0.0277 0.002577 -0.2394
0.3562 -0.1569 -0.6226 0.1604 -0.01731
0.3248 0.432 0.1625 -0.3244 -0.1595
0.3511 0.2785 0.2983 0.2819 0.4118
0.02568 0.3369 -0.374 0.3135 0.4796
0.3167 -0.1565 0.1411 0.4385 0.0182
0.2352 0.2019 -0.1087 0.3636 -0.565
0.3209 -0.4392 0.01931 -0.315 0.4424
0.3052 0.2951 -0.3886 -0.5123 -0.03973
, Matrix(5 x 5)
2.288 1.517 1.607 1.892 1.183
0 1.105 0.7235 0.07972 0.07877
0 0 0.6674 0.299 -0.4158
0 0 0 0.4826 0.6031
0 0 0 0 0.9661
]
is a valid decomposition of the data.
The modified Gram-Schmidt QR decomposition introduced to overcome numerical problems of the Gram-Schmidt QR is given by the algorithm below:
for \(i = 1, \ldots, n\):
\(v_{i} = a_{i}\)
for \(i = 1, \ldots, n\):
\(r_{ii} = \lVert v_{i}\rVert\)
\(q_{i} = v_{i}/r_{ii} \)
for \(j = i + 1, \ldots, n\):
\(r_{ij} = q_{i}^{T}v_{j}\)
\(v_{j} = v_{j} - r_{ij}q_{i}\)
the code for the modified Gram-Schmidt algorithm is given below:
auto qr(QRAlgorithm algo, T)(Matrix!(T) A)
if(isFloatingPoint!(T) && (algo == ModifiedGramSchmidt))
{
auto Q = Matrix!(T)(A.nrow, A.ncol);
auto V = A.dup; long n = A.ncol;
auto R = zeros!(T)(A.ncol, A.ncol);
for(long i = 0; i < n; ++i)
{
T rii = norm(V[i]);
R[i, i] = rii;
Q[i] = map!((a) => a/rii)(V[i]);
for(long j = i + 1; j < n; ++j)
{
T rij = dot(Q[i], V[j]);
R[i, j] = rij;
V[j] = map!((a, b) => a - b*rij)(V[j], Q[i]);
}
}
return [Q, R];
}
run with:
writeln("qr!(ModifiedGramSchmidt)(X): \n", qr!(ModifiedGramSchmidt)(X));
gives the output:
qr!(ModifiedGramSchmidt)(X):
[Matrix(10 x 5)
0.3757 0.1337 0.4163 -0.07128 -0.02322
0.3884 -0.4844 0.0277 0.002577 -0.2394
0.3562 -0.1569 -0.6226 0.1604 -0.01731
0.3248 0.432 0.1625 -0.3244 -0.1595
0.3511 0.2785 0.2983 0.2819 0.4118
0.02568 0.3369 -0.374 0.3135 0.4796
0.3167 -0.1565 0.1411 0.4385 0.0182
0.2352 0.2019 -0.1087 0.3636 -0.565
0.3209 -0.4392 0.01931 -0.315 0.4424
0.3052 0.2951 -0.3886 -0.5123 -0.03973
, Matrix(5 x 5)
2.288 1.517 1.607 1.892 1.183
0 1.105 0.7235 0.07972 0.07877
0 0 0.6674 0.299 -0.4158
0 0 0 0.4826 0.6031
0 0 0 0 0.9661
]
The Householder QR decomposition algorithm is given below:
for \(i = 1, \ldots, n\):
\(x = A_{i:m, i}\)
\(v_{i} = sign(x_{1})\lVert x\rVert e_{1} + x\)
\(v_{i} = v_{i}/\lVert v_{i}\rVert\)
\(A_{i:m, i:n} = A_{i:m, i:n} - 2v_{i}(v_{i}^{T}A_{i:m, i:n})\)
In the above expression, \(x_{1}\) is the first element of the \(x\) vector and \(e_{1}\) is the unit vector for the first element \(\{1, 0, \ldots, 0\}\). The \(Q_{i}\) matrix (at iteration \(i\)) is given by:
$$
Q_{i} = \begin{bmatrix}
I & 0\\
0 & F
\end{bmatrix}
$$
\(I\) is a \((i - 1) \times (i - 1)\) identity matrix and \(F\) is an \((m - i + 1) \times (m - i + 1)\) unitary matrix computable by:
$$ F = I - 2\frac{v v^{T}}{v^{T}v} $$
additionally:
$$ Q = Q_{1}Q_{2} \ldots Q_{n} $$
in this case the returned \(Q\) is an \(m \times m\) matrix and can be computed while \(R\) is being calculated, as in the implementation below:
auto qr(QRAlgorithm algo, T)(Matrix!(T) A)
if(isFloatingPoint!(T) && (algo == Householder))
{
auto _R = A.dup;
long n = A.ncol;
Matrix!(T) Q;
for(long i = 0; i < n; ++i)
{
T[] v = _R[i..$, i].array;
v[0] += norm(v)*sign(v[0]);
T normV = norm(v);
v = map!((a) => a/normV)(v);//later we turn these into Q_i
if(i == 0)
{
Q = F(v);
}else{
updateQ(Q, F(v), i);
}
T[] twoV = map!((a) => 2*a)(v);
_R[i..$, i..$] -= outer(twoV, gemv!(Column)(_R[i..$, i..$], v));
}
auto R = zeros!(T)(n, n);
for(long i = 0; i < n; ++i)
{
for(long j = i; j < n; ++j)
{
R[i, j] = _R[i, j];
}
}
return [Q, R];
}
The functions for calculating \(F\) and updating \(Q\) are given below:
Matrix!(T) F(T)(T[] v)
{
auto twoV = map!((a) => 2*a)(v);
return eye!(T)(v.length) - outer(twoV, v)/norm!(false)(v);
}
auto updateQ(T)(ref Matrix!(T) Q, Matrix!(T) f, size_t i)
{
Matrix!(T) Qk = eye!(T)(i + f.ncol);
Qk[i..$, i..$] = f;
Q = gemm(Q, Qk);
}
In the updateQ() function we call gemm for simplicity but it is computationally more efficient to do the appropriate block-wise matrix multiplication to update Q.
Applying the Householder QR algorithm to our \(X\) matrix:
writeln("qr!(Householder)(X): \n", qr!(Householder)(X));
gives the output:
qr!(Householder)(X):
[Matrix(10 x 10)
-0.3757 0.1337 0.4163 -0.07128 0.02322 0.1172 -0.4304 -0.6402 0.2301 -0.01684
-0.3884 -0.4844 0.0277 0.002577 0.2394 0.4451 -0.178 0.1594 -0.5327 0.1322
-0.3562 -0.1569 -0.6226 0.1604 0.01731 -0.3945 -0.1072 -0.2716 -0.1028 -0.4281
-0.3248 0.432 0.1625 -0.3244 0.1595 0.1787 0.3727 0.1742 -0.1429 -0.5737
-0.3511 0.2785 0.2983 0.2819 -0.4118 -0.4271 -0.193 0.3451 -0.3185 0.1447
-0.02568 0.3369 -0.374 0.3135 -0.4796 0.6392 -0.08517 -0.02471 0.0293 -0.01923
-0.3167 -0.1565 0.1411 0.4385 -0.0182 0.01226 0.7309 -0.2571 0.118 0.2196
-0.2352 0.2019 -0.1087 0.3636 0.565 0.03519 -0.2041 0.4025 0.4777 0.08742
-0.3209 -0.4392 0.01931 -0.315 -0.4424 0.03764 0.002359 0.3289 0.536 -0.1088
-0.3052 0.2951 -0.3886 -0.5123 0.03973 -0.08232 0.1103 -0.06733 -0.02333 0.6171
, Matrix(5 x 5)
-2.288 -1.517 -1.607 -1.892 -1.183
0 1.105 0.7235 0.07972 0.07877
0 0 0.6674 0.299 -0.4158
0 0 0 0.4826 0.6031
0 0 0 0 -0.9661
Which returns the full \(Q\) matrix and reduced \(R\) QR decomposition.
That’s all. Thanks for reading.