Householder transformation
NOTES ON STATISTICS, PROBABILITY and MATHEMATICS
Householder Transformation:
From here:
The more common approach to QR decomposition is employing Householder reflections rather than utilizing Gram-Schmidt. In practice, the Gram-Schmidt procedure is not recommended as it can lead to cancellation that causes inaccuracy of the computation of \(q_j,\) which may result in a non-orthogonal \(Q\) matrix.
Mathematical Concepts:
The projection of \(\vec x\) on \(\vec u\) (unit vector perpendicular to the plane) is a scalar multiple of \(\vec u\), we find that \(\vec w = \vec u \, \alpha\). It follows that \(0 = \vec u^T\,(\vec x - \vec w) = \vec u^T\,(\vec x - \vec u\,\alpha)\). Therefore,
\[\alpha = \frac{\vec u^T \vec x}{\vec u^T \vec u}\]
And substituting,
\[\vec w = \vec u \, \alpha = \frac{\vec u\,\vec u^T}{\vec u^T \vec u} \vec x\]
Since \(\vec u\) is of unit length,
\[\vec w = \vec u\,\vec u^T \, \vec x = \vec u^T \vec x \vec u\]
as on the diagram above. And the reflection of \(\vec x\) is:
\[\vec y = \vec x - 2\,\vec w\\ = \vec x -2\,\vec u^T \vec x \, \vec u \\= \vec x -2\, \vec u \vec u^T \vec x \\= (I - 2\, \vec u \vec u^T)\,\vec x\]
with the matrix \(H = I - 2\,\vec u \vec u^T\) being the mirror or reflection matrix (aka a Householder transformation).
A YT presentation on this topic can be found here.
Since we want to get an upper-triangular matrix \(R\) in the intended \(QR\) decomposition, the computation is carried out by successive sub-matrices such as in each one, the first column is of the form \(\begin{bmatrix}a_1,0,0,\cdots\end{bmatrix}^\top\). Since the \(H\) matrix is orthonormal, the transformation exerted by \(H\vec x\) leaves the norm of \(\vec x\) unchanged. We want a “mirror” such that the projection is along \(e_1\) and of length \(||\vec x||\). Given the equal norm of both \(\vec x\) and \(\vec y\), and with only one non-zero entry in \(\vec y\), \(\vec y= H\vec x = \beta \vec e_1 =\begin{bmatrix} \pm \lVert x\rVert&0&0&\cdots&0\end{bmatrix}\), imagining the reflection along the \(x\) axis results in a simple expression for the difference vector \(\vec v\):
\[\vec v = \vec x - \vec y= \begin{bmatrix}x_1\\x_2\\x_3\\\vdots\\x_m\end{bmatrix} \pm \lVert x \rVert \begin{bmatrix}1\\0\\0\\\vdots\\0\end{bmatrix}\] Interestingly, this norm of \(\vec x\) can be multiplied by \(\pm 1\) or any complex number of norm \(1,\) because it will only change the plane of symmetry of the reflection.
In the actual calculation of \(R\), the procedure involves progressively nested sub-matrices, explaining why there is only one non-zero element in the first column matrix. This can be seen in the code below.
Having figured out \(\vec v\) thus, the unit vector perpendicular to the plane of symmetry will be \(\vec u= \frac{v}{\lVert v \rVert}\).
In the diagram below, the reflected vector is aligned with the \(x\)-axis, and we are trying to figure out the plane of reflection \(\vec u.\)
Source in here
From the diagrams above,
\[\vec x - 2 \vec w^\top \, \vec x \, \vec w = \vec x - 2 \, \vec w \, \vec w \, \vec x = (\color{red}{1 - 2 \, \vec w \, \vec w^\top})\, \vec x\]
In matrix form
\[(\color{red}{I - 2 \, \vec w\, \vec w^\top})\, \vec x\]
In red is the matrix that captures the reflection like on a mirror, the so-called Household transformation or reflector.
What if we are interested in the mirror or dividing hyperplane given two vectors of equal magnitude?
\[\vec x = \vec v + \vec y\]
We want a sequence of unitary matrices such that
\[H_{n-1}\,\dots,H_0\, A = R\]
with \(R\) being an upper triangular matrix:
\[H_0\,\begin{bmatrix}\circ& \cdots & \circ \\ \vdots & \ddots & \circ \\ \circ & \cdots & \circ\end{bmatrix} = \begin{bmatrix}X& \cdots & \circ \\ 0 & \ddots & \circ \\ 0 & \cdots & \circ\end{bmatrix}\]
The matrix \(H_0\) transforms the first column vector into a multiple of the first standard basis \(\vec e_1 = \begin{bmatrix}X & 0 & 0 & \cdots \end{bmatrix}^\top\).
Calling the fist vector \(x,\)
\[H\,x = \beta \, e_0\]
\(H\) preserves length, and therefore, \(\vert \beta \vert = \Vert x\Vert_2.\)
As stated above and apparent on the diagrams:
\[\vec v= \vec x - \vec y = \begin{bmatrix}\chi_0\\\chi_1\\\vdots\\\chi_{n-1}\end{bmatrix}- \begin{bmatrix} \Vert x \Vert \\0\\\vdots\\0\end{bmatrix}\]
Yet, \(\beta\) can be the 2-norm of \(x\), its negative or multiplied by any complex number of norm \(1\).
So we have
\[\vec y =(I - 2 \vec w \vec w^\top)\, \vec x = \beta \,\vec e_0\]
or
\[H\vec x = \pm \Vert x\Vert_2 \, \vec e_0\] Since as above
\[\vec v = x \pm \Vert x\Vert_2 \, \vec e_0 = \begin{bmatrix}\chi_0 \pm \Vert x \Vert_2\\ \chi_1\\ \vdots \\ \chi_{m-1} \end{bmatrix}\]
it’s important to avoid catastrophic cancellation if the weight of the vector \(x\) is concentrated on \(\chi_0\) and we choose \(-\). Therefore, we should make sure that we multiply \(\Vert x\Vert_2\) by the \(\text{sgn} (\chi_0)\).
Instead of the normalizing the vector \(\vec v\) to get \(\vec u = \frac{\vec v}{\Vert v\Vert},\) we can let \(\vec v\) scale, and later normalize \(\vec u\) as in
\[I - \tau\, u\,u^\top\]
with \[\tau = \frac 2{1 + u^\top\, u}\]
We want to scale \(\vec v\) in a way
that \(\chi_0 \pm \Vert x\Vert_2 =1.\)
Now \(\vec u\) doesn’t have to be
unitary. The vector \(\vec u\) is the
Householder vector. Modified (see code below) it is
stored in the lower triangular matrix in the compact output of
qr().
Tau will also be stored in the output qraux as in the
code below.
From the first column of the matrix \(A,\) the Householder vector \(\vec u\) is computed so that
\[(I - \tau_0\, u_0u_0^\top)\, A = \begin{bmatrix}X& \cdots & \circ \\ 0 & \ddots & \circ \\ 0 & \cdots & \circ\end{bmatrix}\]
Next we leave the first row alone, and do the same procedure on the second column of the matrix resulting from multiplying \(H_0A\):
\[\begin{bmatrix}X& \triangle &\cdots & \circ \\ 0 & \square & \cdots & \circ\\0 &\color{red}\triangle & \ddots & \circ \\ 0 & \color{red}\triangle & \cdots & \circ\end{bmatrix} =\begin{bmatrix}X& \triangle &\cdots & \circ \\ 0 & \color{blue}\square & \cdots & 0\\ 0 &\color{red}0 & \ddots & \circ \\ 0 & \color{red}0 & \cdots & \circ\end{bmatrix}\]
The function
\[\text{Houser} \begin{pmatrix}\chi_1 \\ x_2\end{pmatrix} = \begin{bmatrix}\begin{pmatrix} e\\ u_2 \end{pmatrix}, \tau \end{bmatrix}\]
with \(\vec x = \begin{pmatrix}\chi_1 \\ x_2\end{pmatrix}\), \(\begin{pmatrix} e\\ u_2 \end{pmatrix}\) the Householder vector with a \(1\) as the first element. The computation of \(e, u_2, \tau\) is so that
\[I - \tau \, \begin{pmatrix}1 \\ u_2 \end{pmatrix}\begin{pmatrix}1 \\ u_2 \end{pmatrix}^\top \begin{pmatrix}\chi_1 \\ x_2\end{pmatrix} = \begin{pmatrix} e \\ 0\end{pmatrix}\]
Expressed as block matrices:
\[\begin{pmatrix} \bf I & \Bigg \vert & \bf 0 \\ \hline \bf 0 & \Bigg \vert & I - \tau_1 \, \begin{pmatrix}1\\u_{21}\end{pmatrix}\begin{pmatrix}1\\u_{21}\end{pmatrix}^\top \end{pmatrix} \, \begin{pmatrix} \bf R_{00} & \Bigg \vert & r_{01} & R_{02} \\ \hline \\ & \Bigg \vert & \alpha_{11} & a_{12}^\top\\ \bf 0 & \Bigg \vert \\ & \Bigg \vert & a_{21} & A_{22} \end{pmatrix}\]
Computation:
The function House() below has embedded some
print() lines to reproduce the example elaborated in here.
House <- function(A) {
require(Matrix)
R <- as.matrix(A) # Set R to the input matrix A
n <- ncol(A)
m <- nrow(A)
H <- list() # Initialize a list to store the computed H matrices to calculate Q later
if (m > n) {c <- n} else {c <- m}
for (k in 1:c) {
x <- R[k:m, k] # Equivalent to a_1
print(paste0("The HR for iteration ", k)); print(round(R, 2))
print(paste0("The x for iteration ", k)); print(x)
print(paste0("The x1 of the column matrix is ", x[1], " sign ", sign(x[1])))
e1 <- as.matrix(c(1, rep(0, length(x) - 1))) # e_1
vk <- x + sign(x[1]) * sqrt(sum(x^2)) * e1 # v as x + y
print("The reflecion vec v = x + norm x e_1 is: "); print(x); print(" + "); print(as.vector(sign(x[1]) * sqrt(sum(x^2)) * e1)); print(" = "); print(as.vector(vk))
# Compute the H matrix
hk <- diag(length(x)) - 2 * as.vector(vk %*% t(vk)) / (t(vk) %*% vk) # I - 2 vv^T / vTv
print("vv^T: "); print(round(vk %*% t(vk), 2))
print("The hk is:"); print(round(hk, 2))
if (k > 1) {# After the first matrix H a diagonal of 1's is added
hk <- bdiag(diag(k - 1), hk)
print("The hk matrix with 1's is:")
print(hk)
}
# Store the H matrix to find Q at the end of iteration
H[[k]] <- hk
R <- hk %*% R # The first R is the actual A and it is updated as
# H H H H A
}
Q <- Reduce("%*%", H) # Calculate Q matrix by multiplying all H matrices
res <- list('Q'= round(Q, 2),'R'= round(R, 2))
return(res)
}Now with the input from here:
\[A = \begin{bmatrix} 2 & – 2 & 18 \\\ 2 & 1 & 0 \\\ 1 & 2 & 0 \end{bmatrix}\]
The first column vector, i.e. \(a_1 = \begin{bmatrix}2 & 2 & 1 \end{bmatrix}^T\) will be reflected to \(\vec x - \vec y\). But we can instead of using \(-||x||\) we can use that value times \(-1\) and calculate \(\vec x + ||x||\,\vec e_1\):
\[\vec v_1 = \begin{bmatrix}2\\2\\1\end{bmatrix} + \left( \sum_{k=1}^m a_1^2 \right)^{1/2} \, \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix} = \begin{bmatrix}5 \\ 2 \\ 1\end{bmatrix}\]
The next step is calculating \(I - 2 \frac{vv^\top}{v^\top v}\):
\[H_1= \small \begin{bmatrix}1&0&0\\0&1&0\\0&0&1\end{bmatrix} - 2 \frac{\begin{bmatrix}5&2&1\end{bmatrix}\begin{bmatrix}5\\2\\1\end{bmatrix}}{\begin{bmatrix}5\\2\\1\end{bmatrix}\begin{bmatrix}5&2&1\end{bmatrix}}=\begin{bmatrix} -\frac{2}{3} & -\frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & 0.7333 & -0.1333 \\\ -\frac{1}{3} & -0.1333 & 0.9333 \end{bmatrix}\]
This is then multiplied \(H_1A\):
\[H_1 A = \small \begin{bmatrix} -\frac{2}{3} & -\frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & 0.7333 & -0.1333 \\\ -\frac{1}{3} & -0.1333 & 0.9333 \end{bmatrix} \begin{bmatrix} 2 & – 2 & 18 \\\ 2 & 1 & 0 \\\ 1 & 2 & 0 \end{bmatrix} = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & 1.8 & 12 \\\ 0 & 2.4 & -6 \end{bmatrix}\]
\(H_1A\) replaces \(A\) in the next iteration. Now comes to \(a_2\) and \(H_2\) with the submatrix of \(A\)
\[A^{(1)} = \begin{bmatrix} 1.8 & 12 \\\ 2.4 & -6 \end{bmatrix}\] \(v_2\) is
\[\small \begin{bmatrix} 1.8 \\\ 2.4 \end{bmatrix} + \sqrt{\sum^m_{j=1} a_2^2} \begin{bmatrix} 1 \\\ 0 \end{bmatrix} = \begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix}\]
The \(H_2\) Householder matrix is
\[H_2 = \small \begin{bmatrix} 1 & 0 \\\ 0 & 1 \end{bmatrix} – 2 \frac{\begin{bmatrix} 4.8 & 2.4 \end{bmatrix} \begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix}}{\begin{bmatrix} 4.8 \\\ 2.4 \end{bmatrix} \begin{bmatrix} 4.8 & 2.4 \end{bmatrix}}=\begin{bmatrix} 1 & 0 & 0 \\\ 0 & -0.6 & -0.8 \\\ 0 & -0.8 & 0.6 \end{bmatrix}\]
The first column \(\begin{bmatrix}1 \\\ 0 \\\ 0 \end{bmatrix}\) and first row \(\begin{bmatrix}1 & 0 & 0 \end{bmatrix}\) are added to the resulting \(H_2\) matrix to keep it \(n \times n\).
Now
\[H_2 = \begin{bmatrix} 1 & 0 & 0 \\\ 0 & -0.6 & -0.8 \\\ 0 & -0.8 & 0.6 \end{bmatrix} \begin{bmatrix} -3 & 0 & -12 \\\ 0 & 1.8 & 12 \\\ 0 & 2.4 & -6 \end{bmatrix}=\begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & 6 \end{bmatrix}\]
Moving to \(a_3\) and \(H_3\), the submatrix of \(H_2A\) is simply \([6]\). Therefore, \(v_3\) is equal to:
\[\begin{bmatrix} 6 \end{bmatrix} – \small \sqrt{\sum^m_{j=1} a_3^2} \begin{bmatrix} 1 \end{bmatrix} = 12\]
The \(H_3\)
\[H_3 = \small \begin{bmatrix} 1 \end{bmatrix} – 2 \frac{\begin{bmatrix} 12 \end{bmatrix}\begin{bmatrix} 12 \end{bmatrix}}{\begin{bmatrix} 12 \end{bmatrix}\begin{bmatrix} 12 \end{bmatrix}} = \begin{bmatrix}1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & -1 \end{bmatrix}\]
and
\[H_3 A = \small \begin{bmatrix}1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & -1 \end{bmatrix} \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & 6 \end{bmatrix} = \begin{bmatrix} -3 & 0 & -12 \\\ 0 & -3 & 12 \\\ 0 & 0 & -6 \end{bmatrix}\]
Which is the \(R\) factorization in the \(QR\) decomposition method. The \(Q\) factorization of \(QR\) decomposition is found by multiplying all the \(H\) matrices together
\[Q = \small \begin{bmatrix} -\frac{2}{3} & -\frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & 0.7333 & -0.1333 \\\ -\frac{1}{3} & -0.1333 & 0.9333 \end{bmatrix} \begin{bmatrix} 1 & 0 & 0 \\\ 0 & -0.6 & -0.8 \\\ 0 & -0.8 & 0.6 \end{bmatrix} \begin{bmatrix}1 & 0 & 0 \\\ 0 & 1 & 0 \\\ 0 & 0 & -1 \end{bmatrix}= \begin{bmatrix} -\frac{2}{3} & \frac{2}{3} & -\frac{1}{3} \\\ -\frac{2}{3} & -\frac{1}{3} & \frac{2}{3} \\\ -\frac{1}{3} & -\frac{2}{3} & -\frac{2}{3} \end{bmatrix}\]
A <- rbind(c(2,-2,18),c(2,1,0),c(1,2,0))
(hse <- House(A))## [1] "The HR for iteration 1"
## [,1] [,2] [,3]
## [1,] 2 -2 18
## [2,] 2 1 0
## [3,] 1 2 0
## [1] "The x for iteration 1"
## [1] 2 2 1
## [1] "The x1 of the column matrix is 2 sign 1"
## [1] "The reflecion vec v = x + norm x e_1 is: "
## [1] 2 2 1
## [1] " + "
## [1] 3 0 0
## [1] " = "
## [1] 5 2 1
## [1] "vv^T: "
## [,1] [,2] [,3]
## [1,] 25 10 5
## [2,] 10 4 2
## [3,] 5 2 1
## [1] "The hk is:"
## [,1] [,2] [,3]
## [1,] -0.67 -0.67 -0.33
## [2,] -0.67 0.73 -0.13
## [3,] -0.33 -0.13 0.93
## [1] "The HR for iteration 2"
## [,1] [,2] [,3]
## [1,] -3 0.0 -12
## [2,] 0 1.8 -12
## [3,] 0 2.4 -6
## [1] "The x for iteration 2"
## [1] 1.8 2.4
## [1] "The x1 of the column matrix is 1.8 sign 1"
## [1] "The reflecion vec v = x + norm x e_1 is: "
## [1] 1.8 2.4
## [1] " + "
## [1] 3 0
## [1] " = "
## [1] 4.8 2.4
## [1] "vv^T: "
## [,1] [,2]
## [1,] 23.04 11.52
## [2,] 11.52 5.76
## [1] "The hk is:"
## [,1] [,2]
## [1,] -0.6 -0.8
## [2,] -0.8 0.6
## [1] "The hk matrix with 1's is:"
## 3 x 3 sparse Matrix of class "dsCMatrix"
##
## [1,] 1 . .
## [2,] . -0.6 -0.8
## [3,] . -0.8 0.6
## [1] "The HR for iteration 3"
## 3 x 3 Matrix of class "dgeMatrix"
## [,1] [,2] [,3]
## [1,] -3 0 -12
## [2,] 0 -3 12
## [3,] 0 0 6
## [1] "The x for iteration 3"
## [1] 6
## [1] "The x1 of the column matrix is 6 sign 1"
## [1] "The reflecion vec v = x + norm x e_1 is: "
## [1] 6
## [1] " + "
## [1] 6
## [1] " = "
## [1] 12
## [1] "vv^T: "
## [,1]
## [1,] 144
## [1] "The hk is:"
## [,1]
## [1,] -1
## [1] "The hk matrix with 1's is:"
## 3 x 3 sparse Matrix of class "dsCMatrix"
##
## [1,] 1 . .
## [2,] . 1 .
## [3,] . . -1## $Q
## 3 x 3 Matrix of class "dgeMatrix"
## [,1] [,2] [,3]
## [1,] -0.67 0.67 -0.33
## [2,] -0.67 -0.33 0.67
## [3,] -0.33 -0.67 -0.67
##
## $R
## 3 x 3 Matrix of class "dgeMatrix"
## [,1] [,2] [,3]
## [1,] -3 0 -12
## [2,] 0 -3 12
## [3,] 0 0 -6Now comparing to built-in functions:
round(qr.Q(qr(A)), 2); round(qr.R(qr(A)), 2)## [,1] [,2] [,3]
## [1,] -0.67 0.67 0.33
## [2,] -0.67 -0.33 -0.67
## [3,] -0.33 -0.67 0.67## [,1] [,2] [,3]
## [1,] -3 0 -12
## [2,] 0 -3 12
## [3,] 0 0 6Analyzing the result of qr(), note that qr
is a matrix with the same dimensions as the input x. The
upper triangle contains the \(R\) of
the decomposition and the lower triangle contains information on the
\(Q\) of the decomposition.
(QR <- round(as.matrix(qr(A)$qr), 2))## [,1] [,2] [,3]
## [1,] -3.00 0.0 -12
## [2,] 0.67 -3.0 12
## [3,] 0.33 0.8 6# Excluding the lower triangular entries returns R:
QR[row(QR) > col(QR)] = 0; QR## [,1] [,2] [,3]
## [1,] -3 0 -12
## [2,] 0 -3 12
## [3,] 0 0 6# To get Q:
round(Q <- qr.Q(qr(A)), 2)## [,1] [,2] [,3]
## [1,] -0.67 0.67 0.33
## [2,] -0.67 -0.33 -0.67
## [3,] -0.33 -0.67 0.67# Notice that Q is orthonormal:
round(t(Q)%*%Q, 2)## [,1] [,2] [,3]
## [1,] 1 0 0
## [2,] 0 1 0
## [3,] 0 0 1To replicate the lower triangular matrix in the output of R for
qr() as well as qraux:
Source in here
Formulas:
qraux entries contain the tau values
obtained for each H sub-matrix:
\[\tau = \frac{2}{1 + \vec u^\top\,\vec u}\]
such that it is the scalar factor in the Householder matrix:
\[I - \tau \,\vec u\, \vec u^\top\]
and the vector \(\vec u\) is the Householder vector from the leading column:
\[\vec u = x /\nu\]
with
\[\nu = \chi_1 + \text{sign}(\chi_1)\,\Vert x\Vert_2 = \chi_1 + \text{sign}(\chi_1)\, \sqrt{\sum_i x_i^2}\]
where \(\chi_1\) is the leading entry in the column vector, and \(\Vert x \Vert_2\) the square root of the 2-norm of the leading column vector in the sub-matrices.
set.seed(2024)
options(digits=10)
# Since it's about understanding the convention of the qr()
# I'm sticking to square matrices. s is the number of rows and cols:
s = 4
m = matrix(rnorm(s^2), nrow = s)
# chi1 is the leading entry of each col vector in the matrix m (aka as A)
# keeping in mind that the matrix will be updated in a loop
# by multiplying it on the left by the Householder matrices:
chi1 = m[1,1] # This is the initial entry out of the original matrix.
# The initial col vector is precisely the first col in the orig. mat.
# And its 2-norm is:
norm = sqrt(sum(m[,1]^2))
# The Householder vector will be constructed by taking the entries of
# the column matrix at each iteration of the products H H H A
# below the leading entry in the decreasing sub-matrices (chi1):
x2 = m[2:nrow(m), 1] # This is the initial one.
# nu1 is the difference between the original vector and the mirror image
# avoiding catastrophic cancellation, ie chi1 + norm(x) * e_1:
nu1 = chi1 + sign(chi1) * norm
# nu1 will standardize u2:
u2 = x2 / nu1
# Let's see this applied to this first column vector in the matrix A
(tau = 2 / (1 + t(u2)%*%u2))## [,1]
## [1,] 1.881503125# Now let's get all the values of the output in qraux and in the lower
# triangular matrix of the output:
# Initializing a vector to collect values below the diag
# They will be output with a sign difference wrt to the R base function.
lwrT = matrix(0, s, s)
# Initializing an empty vector to collect the values of qraux:
qraux <- numeric()
A = m
r = nrow(m)
for (i in 1:r){
chi1 = A[i,i]
norm = sqrt(sum(A[i:r, i]^2))
if(i < r){x2 = as.matrix(A[(i + 1):r, i])}else{x2 = 0}
nu1 = chi1 + sign(chi1) * norm
u2 = x2 / nu1
if(i < r){tau = 2 / (1 + t(u2) %*% u2)}else{tau = sign(chi1) * nu1 / 2}
qraux[i] = tau
if(i < r) {augm = c(rep(0, i - 1), 1, u2)}
else
{augm = c(rep(0, i - 1), 0)}
H <- diag(r) - as.numeric(tau) *
outer(as.vector(augm), as.vector(t(augm)),'*')
if(i < r){lwrT[(i + 1):r, i] <- - H[(i + 1):r, i]}else{}
A <- H %*% A
}
# After all iterations multiplying H A we end up producing
# R or the upper triangular (including the diagonal in qr()):
round(A,9)## [,1] [,2] [,3] [,4]
## [1,] -1.11397156 -1.537094686 1.479110318 -1.202532687
## [2,] 0.00000000 -0.975037638 1.587258014 2.295022409
## [3,] 0.00000000 0.000000000 1.032380986 2.124104691
## [4,] 0.00000000 0.000000000 0.000000000 -0.544290993# The lower triangular containing essentially tau * uTu:
lwrT## [,1] [,2] [,3] [,4]
## [1,] 0.0000000000 0.0000000000 0.0000000000 0
## [2,] 0.4207603289 0.0000000000 0.0000000000 0
## [3,] -0.0969246757 0.6907309238 0.0000000000 0
## [4,] -0.1910983876 0.1504635437 -0.3677919994 0# They fit perfectly to replicate the output of qr():
qr(m)$qr## [,1] [,2] [,3] [,4]
## [1,] -1.1139715602 -1.5370946856 1.4791103176 -1.2025326872
## [2,] 0.4207603289 -0.9750376385 1.5872580139 2.2950224086
## [3,] -0.0969246757 0.6907309238 1.0323809861 2.1241046912
## [4,] -0.1910983876 0.1504635437 -0.3677919994 -0.5442909928# Likewise the output of the vector in the calculations above
# for tau:
qraux## [1] 1.8815031249 1.7072846053 1.9299080843 0.5442909928# replicate the built-in function qraux:
qr(m)$qraux## [1] 1.8815031249 1.7072846053 1.9299080843 0.5442909928NOTE: These are tentative notes on different topics for personal use - expect mistakes and misunderstandings.