Skip to content

Commit 4ba2b77

Browse files
feat(matrix): add QR decomposition algorithm using Gram-Schmidt process
Decomposes a matrix A into an orthogonal matrix Q and an upper triangular matrix R such that A = Q * R. ### What This Adds **QRDecomposition.java** - Main algorithm implementation: - `decompose()` - Performs QR factorization using the Gram-Schmidt process - Returns a QR object containing both Q (orthogonal) and R (upper triangular) matrices - Validates input matrix using MatrixUtil.validateInputMatrix() - Throws ArithmeticException for rank-deficient matrices **QRDecompositionTest.java** - Unit tests: - Tests for 2x2 and 3x3 matrix decomposition - Verifies Q * R reconstruction equals original matrix - Validates Q columns are orthonormal - Confirms R is upper triangular - Tests identity matrix decomposition - Tests rank-deficient matrix rejection - Tests null and empty matrix rejection ### Algorithm The Gram-Schmidt process orthogonalizes columns iteratively: - For each column j, subtract projections onto previous orthogonal vectors - Normalize to get j-th column of Q - Store coefficients in R[i][j] Time: O(m*n^2) | Space: O(m*n + n^2) ### Reference https://en.wikipedia.org/wiki/QR_decomposition
1 parent 4b8099c commit 4ba2b77

2 files changed

Lines changed: 230 additions & 0 deletions

File tree

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
package com.thealgorithms.matrix;
2+
3+
/**
4+
* @brief Implementation of QR Decomposition using the Gram-Schmidt process
5+
* @details Decomposes a matrix A into an orthogonal matrix Q and an upper
6+
* triangular matrix R such that A = Q * R. The Gram-Schmidt process
7+
* orthogonalizes the columns of A to produce Q, and R is computed as Q^T * A.
8+
* This decomposition is useful for solving linear least squares problems,
9+
* eigenvalue computations, and numerical stability in linear algebra.
10+
* @see <a href="https://en.wikipedia.org/wiki/QR_decomposition">QR Decomposition</a>
11+
*/
12+
public final class QRDecomposition {
13+
14+
private QRDecomposition() {
15+
}
16+
17+
/**
18+
* A helper class to store both Q and R matrices
19+
*/
20+
public static class QR {
21+
private final double[][] q;
22+
private final double[][] r;
23+
24+
QR(double[][] q, double[][] r) {
25+
this.q = q;
26+
this.r = r;
27+
}
28+
29+
public double[][] getQ() {
30+
return q;
31+
}
32+
33+
public double[][] getR() {
34+
return r;
35+
}
36+
}
37+
38+
/**
39+
* @brief Performs QR decomposition on a matrix using the Gram-Schmidt process
40+
* @param matrix the input matrix (m x n)
41+
* @return QR object containing orthogonal matrix Q (m x n) and upper triangular matrix R (n x n)
42+
* @throws IllegalArgumentException if the matrix is null, empty, or has invalid rows
43+
*/
44+
public static QR decompose(double[][] matrix) {
45+
MatrixUtil.validateInputMatrix(matrix);
46+
47+
int m = matrix.length;
48+
int n = matrix[0].length;
49+
50+
double[][] q = new double[m][n];
51+
double[][] r = new double[n][n];
52+
53+
for (int j = 0; j < n; j++) {
54+
double[] v = getColumn(matrix, j);
55+
56+
for (int i = 0; i < j; i++) {
57+
double[] qi = getColumn(q, i);
58+
r[i][j] = dotProduct(qi, v);
59+
v = subtractVectors(v, scalarMultiply(qi, r[i][j]));
60+
}
61+
62+
r[j][j] = norm(v);
63+
if (r[j][j] == 0) {
64+
throw new ArithmeticException("Matrix is rank deficient. Cannot perform QR decomposition.");
65+
}
66+
double[] qj = scalarMultiply(v, 1.0 / r[j][j]);
67+
setColumn(q, j, qj);
68+
}
69+
70+
return new QR(q, r);
71+
}
72+
73+
private static double[] getColumn(double[][] matrix, int col) {
74+
int m = matrix.length;
75+
double[] column = new double[m];
76+
for (int i = 0; i < m; i++) {
77+
column[i] = matrix[i][col];
78+
}
79+
return column;
80+
}
81+
82+
private static void setColumn(double[][] matrix, int col, double[] column) {
83+
for (int i = 0; i < matrix.length; i++) {
84+
matrix[i][col] = column[i];
85+
}
86+
}
87+
88+
private static double dotProduct(double[] a, double[] b) {
89+
double sum = 0;
90+
for (int i = 0; i < a.length; i++) {
91+
sum += a[i] * b[i];
92+
}
93+
return sum;
94+
}
95+
96+
private static double[] subtractVectors(double[] a, double[] b) {
97+
double[] result = new double[a.length];
98+
for (int i = 0; i < a.length; i++) {
99+
result[i] = a[i] - b[i];
100+
}
101+
return result;
102+
}
103+
104+
private static double[] scalarMultiply(double[] v, double scalar) {
105+
double[] result = new double[v.length];
106+
for (int i = 0; i < v.length; i++) {
107+
result[i] = v[i] * scalar;
108+
}
109+
return result;
110+
}
111+
112+
private static double norm(double[] v) {
113+
return Math.sqrt(dotProduct(v, v));
114+
}
115+
}
Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
package com.thealgorithms.matrix;
2+
3+
import static org.junit.jupiter.api.Assertions.assertArrayEquals;
4+
import static org.junit.jupiter.api.Assertions.assertThrows;
5+
6+
import org.junit.jupiter.api.Test;
7+
8+
public class QRDecompositionTest {
9+
10+
private static final double DELTA = 1e-9;
11+
12+
@Test
13+
public void testQRDecomposition2x2() {
14+
double[][] matrix = {{12, -51}, {6, 167}};
15+
QRDecomposition.QR qr = QRDecomposition.decompose(matrix);
16+
double[][] q = qr.getQ();
17+
double[][] r = qr.getR();
18+
19+
double[][] reconstructed = multiplyMatrices(q, r);
20+
for (int i = 0; i < matrix.length; i++) {
21+
assertArrayEquals(matrix[i], reconstructed[i], DELTA);
22+
}
23+
}
24+
25+
@Test
26+
public void testQRDecomposition3x3() {
27+
double[][] matrix = {{1, 1, 0}, {1, 0, 1}, {0, 1, 1}};
28+
QRDecomposition.QR qr = QRDecomposition.decompose(matrix);
29+
double[][] q = qr.getQ();
30+
double[][] r = qr.getR();
31+
32+
double[][] reconstructed = multiplyMatrices(q, r);
33+
for (int i = 0; i < matrix.length; i++) {
34+
assertArrayEquals(matrix[i], reconstructed[i], DELTA);
35+
}
36+
}
37+
38+
@Test
39+
public void testQROrthogonalColumns() {
40+
double[][] matrix = {{1, 1, 0}, {1, 0, 1}, {0, 1, 1}};
41+
QRDecomposition.QR qr = QRDecomposition.decompose(matrix);
42+
double[][] q = qr.getQ();
43+
44+
for (int i = 0; i < q[0].length; i++) {
45+
for (int j = i; j < q[0].length; j++) {
46+
double dot = 0;
47+
for (int k = 0; k < q.length; k++) {
48+
dot += q[k][i] * q[k][j];
49+
}
50+
if (i == j) {
51+
assertArrayEquals(new double[] {1.0}, new double[] {dot}, DELTA);
52+
} else {
53+
assertArrayEquals(new double[] {0.0}, new double[] {dot}, DELTA);
54+
}
55+
}
56+
}
57+
}
58+
59+
@Test
60+
public void testRIsUpperTriangular() {
61+
double[][] matrix = {{12, -51}, {6, 167}};
62+
QRDecomposition.QR qr = QRDecomposition.decompose(matrix);
63+
double[][] r = qr.getR();
64+
65+
for (int i = 1; i < r.length; i++) {
66+
for (int j = 0; j < i; j++) {
67+
assertArrayEquals(new double[] {0.0}, new double[] {r[i][j]}, DELTA);
68+
}
69+
}
70+
}
71+
72+
@Test
73+
public void testQRDecompositionIdentityMatrix() {
74+
double[][] matrix = {{1, 0, 0}, {0, 1, 0}, {0, 0, 1}};
75+
QRDecomposition.QR qr = QRDecomposition.decompose(matrix);
76+
double[][] q = qr.getQ();
77+
double[][] r = qr.getR();
78+
79+
for (int i = 0; i < matrix.length; i++) {
80+
assertArrayEquals(matrix[i], q[i], DELTA);
81+
assertArrayEquals(matrix[i], r[i], DELTA);
82+
}
83+
}
84+
85+
@Test
86+
public void testQRDecompositionRankDeficientThrows() {
87+
double[][] matrix = {{1, 2}, {2, 4}};
88+
assertThrows(ArithmeticException.class, () -> QRDecomposition.decompose(matrix));
89+
}
90+
91+
@Test
92+
public void testQRDecompositionNullMatrixThrows() {
93+
assertThrows(IllegalArgumentException.class, () -> QRDecomposition.decompose(null));
94+
}
95+
96+
@Test
97+
public void testQRDecompositionEmptyMatrixThrows() {
98+
assertThrows(IllegalArgumentException.class, () -> QRDecomposition.decompose(new double[0][0]));
99+
}
100+
101+
private static double[][] multiplyMatrices(double[][] a, double[][] b) {
102+
int m = a.length;
103+
int n = b[0].length;
104+
int k = a[0].length;
105+
double[][] result = new double[m][n];
106+
for (int i = 0; i < m; i++) {
107+
for (int j = 0; j < n; j++) {
108+
for (int p = 0; p < k; p++) {
109+
result[i][j] += a[i][p] * b[p][j];
110+
}
111+
}
112+
}
113+
return result;
114+
}
115+
}

0 commit comments

Comments
 (0)