Article directory
- Introduction
- Prepositional theory
- Householder performs QR decomposition
- Code
Introduction
Elementary transformation tools such as trigonometric decomposition (LU decomposition) can be used to solve systems of linear equations, but they do have some limitations. For example, for ill-conditioned systems of linear equations, LU decomposition may lead to numerically unstable results. Furthermore, for irreversible matrices, LU decomposition is not applicable.
To overcome these problems, QR decomposition is introduced, where the matrix is decomposed into an orthogonal matrix Q and an upper triangular matrix R. QR decomposition is applicable to any invertible matrix and can provide numerically stable solutions. The implementation of QR decomposition can be accomplished with the help of techniques such as Schmidt orthogonal normalization, Givens transform and Householder transform.
In QR decomposition, the columns of the orthogonal matrix Q are orthogonal, meaning that they satisfy that the transpose of Q times Q equals the identity matrix. This property helps reduce the propagation of numerical errors and provides numerical stability. At the same time, the upper triangular matrix R contains important information of the original matrix and can be used for tasks such as solving linear equations.
TheHouseholdertransformation is a linear transformation that projects a vector into a new direction, usually on a hyperplane. It can be used to zero out all but the first element in a vector, which allows us to use it for reflection operations in QR decomposition. By successively applying a series of Householder transformations, we can transform the original matrix A into an upper triangular matrix R. In this process, we also construct the orthogonal matrix Q, which will be used for the final QR decomposition.
The key idea ofHouseholder transformation is to find a reflective surface (or reflective hyperplane) that can map a vector to a zero vector or a specific direction. The design of this transformation allows us to zero out certain elements, thus achieving the upper triangular shape of R. Householder transform is widely used in numerical calculations and linear algebra, especially in fields such as QR decomposition and eigenvalue calculation.
Pre-positioned theory
Let’s look at a theorem first:
For any vector with norm 1
ω
∈
R
n
\omega \in {R^n}
ω∈Rn, its reflection matrix is:
H
=
I
?
2
ω
ω
T
H = I – 2\omega {\omega ^T}
H=I?2ωωT, where I is the unit matrix. Reflection array H satisfies:
H
T
=
H
H^{\rm{T}} = H
HT=H
H
?
H
=
I
H*H = I
H?H=I
Brief certificate:
H
T
=
(
I
?
2
ω
ω
T
)
T
=
I
?
2
[
(
ω
T
)
T
ω
T
]
=
I
?
2
ω
ω
T
{H^T} = {\left( {I – 2\omega {\omega ^T}} \right)^T} = I – 2\left[ {{ {\left( {{\omega ^T}} \right)}^T}{\omega ^T}} \right] = I – 2\omega {\omega ^T}
HT=(I?2ωωT)T=I?2[(ωT)TωT]=I?2ωωT
H
?
H
=
(
I
?
2
ω
ω
T
)
?
(
I
?
2
ω
ω
T
)
=
I
?
2
ω
ω
T
?
2
ω
ω
T
+
4
ω
(
ω
T
ω
)
ω
T
=
I
?
4
ω
ω
T
+
4
ω
ω
T
=
I
H * H = \left( {I – 2\omega {\omega ^T}} \right) * \left( {I – 2\omega {\omega ^T}} \right ) = I – 2\omega {\omega ^T} – 2\omega {\omega ^T} + 4\omega \left( {{\omega ^T }\omega } \right){\omega ^T} = I – 4\omega {\omega ^T} + 4\omega {\omega ^T} = I
H?H=(I?2ωωT)?(I?2ωωT)=I?2ωωT?2ωωT + 4ω(ωTω)ωT=I?4ωωT + 4ωωT=I
Reflection transformation can be understood as a process in which two vectors are symmetrical about a normal plane. Assume that the hyperplane S (passes the origin, with w as the normal vector) is:
S
=
{
x
∣
ω
T
x
=
0
,
?
x
∈
R
n
}
S = \{ x{|}{\omega ^T}x = 0,{\forall _x} \in {R^n}\}
S={x∣ωTx=0,?x?∈Rn}
That is:
?
z
,
z
‘
∈
R
n
\forall z，z’ \in {R^n}
?z, z′∈Rn If the two norms of the two vectors are equal, then there is a reflection transformation matrix H such that
z
‘
=
H
z
z’ = Hz
z′=Hz z′ and z are symmetric about the hyperplane S, as shown in the figure below.
This point is important, so let’s express it in a simple way: In an n-dimensional space, as long as the modules of two vectors are equal, you can find a hyperplane that makes the two vectors symmetrical about the hyperplane, that is Can find an H to make the above formula true.
Householder performs QR decomposition
First, we expect to decompose any real matrix A into QR, that is, A = Q*R, where Q is a symmetric orthogonal matrix and R is an upper triangular matrix.
Assume that the matrix A to be decomposed is as follows:
Yes
For each column vector v in A, we can find a unit vector such that
v
→
=
α
e
→
\mathop v\limits^ \to = \alpha \mathop e\limits^ \to
v→=αe→, then we can find H such that
H
v
→
=
α
e
→
H\overrightarrow v \ = \alpha \overrightarrow e
Hv
H
n
.
.
.
.
H
3
H
2
H
1
A
=
R
{H_n}….{H_3}{H_2}{H_1}A = R
Hn?….H3?H2?H1?A=R
Since H is symmetric and orthogonal, then
Q
T
=
H
n
.
.
.
.
H
1
{Q^T} = {H_n}….{H_1}
QT=Hn?….H1?, Q is also symmetrical and orthogonal, then
A
=
Q
R
A = \ Q R
A=QR
Code implementation
Algorithm implementation process: (The Kth transformation matrix H is solved as follows)
v
K
→
=
A
K
(
A
K
is a subarray with a length of
(
n
?
k
+
1
)
)
v
K
→
=
v
K
→
?
∥
v
k
→
∥
?
e
k
→
v
k
→
=
v
k
→
∥
v
∥
H
k
=
I
?
2
v
k
?
v
k
T
\overrightarrow {{\ v _K}} = {A_K}\left( {{A_K} is a subarray with length (n – k + 1) } \right)\\overrightarrow {{\ v _K}} = \overrightarrow {{\ v _K}} – \left \| {\overrightarrow {{\ v _k}} } \right\| \cdot \overrightarrow {{e_k}}\\ \\overrightarrow {{\ v _k}} = \frac{{\overrightarrow {{\ v _k}} }}{{\left\| \ v \right\|}}\\\ {H_k} = I – 2{\ v _k} * {\ \v _k}^T
vK?
#include <iostream>
#include <vector>
#include <Eigen/Dense>
#include <Eigen/Core>
void householderQR(Eigen::MatrixXd & amp;A, Eigen::MatrixXd & amp;Q, Eigen::MatrixXd & amp;R) {<!-- -->
int m = A.rows();
int n = A.cols();
Q = Eigen::MatrixXd::Identity(m, m);
R = A;
for (int k = 0; k < n; + + k) {<!-- -->
Eigen::VectorXd x = R.block(k, k, m - k, 1);
Eigen::VectorXd e1 = Eigen::VectorXd::Zero(m - k);
e1(0) = 1;
Eigen::VectorXd v = x - x.norm() * e1;
v /= v.norm();
Eigen::MatrixXd H = Eigen::MatrixXd::Identity(m, m);
H.block(k, k, m - k, m - k) -= 2.0 * v * v.transpose();
R = H * R;
Q = Q * H.transpose();
}
}
int main() {<!-- -->
Eigen::MatrixXd A(3, 3); // Create a 3x3 example matrix
//Fill matrix A
A << 12, -51, 4,
6, 167, -68,
-4, 24, -41;
Eigen::MatrixXd Q, R;
householderQR(A, Q, R);
std::cout << "Matrix Q:\\
" << Q << "\\
\\
";
std::cout << "Matrix R:\\
" << R << "\\
\\
";
Eigen::MatrixXd reconstructed_A = Q * R;
std::cout << "Reconstructed matrix A:\\
" << reconstructed_A << "\\
\\
";
return 0;
}