6.5.5 QR Decomposition and Least Squares

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188

/*

The QR decomposition factors an $m$ by $n$ matrix $A$ with $m \geq n$ into $A = QR$, where $Q$ has
orthonormal columns and $R$ is upper triangular. Orthonormal columns are what make the factorization
useful: $Q$ preserves lengths, so multiplying by $Q^T$ rewrites a problem in a new basis without
distorting distances, and a triangular system is solved by back substitution.

This implementation uses Householder reflections. A Householder reflection is the mirror across the
hyperplane orthogonal to a chosen vector $v$, written $I - 2vv^T/(v^Tv)$, and $v$ is chosen so that
the reflection maps the part of a column below the diagonal onto the axis, zeroing it in one step.
Applying one reflection per column leaves $R$, and accumulating the reflections gives $Q$.
Reflecting is stable because it never scales anything: the alternative of orthogonalizing column by
column with Gram-Schmidt loses orthogonality as the columns approach linear dependence, whereas
reflections keep it to within rounding error.

The overdetermined system $Ax = b$ usually has no solution, so least squares asks for the $x$ whose
residual is shortest. Because $Q$ preserves lengths, $\|Ax - b\|$ equals $\|Rx - Q^Tb\|$, and the
entries of $Q^Tb$ below row $n$ cannot be changed by any $x$; minimizing therefore means solving the
square triangular system formed by the first $n$ rows, which back substitution does directly.

- `qr_decompose(a, &q, &r)` factors the $m$ by $n$ matrix `a` with $m \geq n$ into an $m$ by $m$
  orthogonal matrix `q` and an $m$ by $n$ upper triangular matrix `r`, whose product is `a`.
- `least_squares(a, b, eps = 1e-10)` returns the vector $x$ minimizing the Euclidean norm of
  $`a`x - `b`$, or `std::nullopt` if the columns of `a` are linearly dependent to within `eps`. The
  matrix `a` must have at least as many rows as columns, and `b` must have one entry per row.

Fitting any linear model reduces to this: one row per data point, one column per basis function.
Solving the normal equations $A^TAx = A^Tb$ by the row reduction of section 6.5.2 gives the same
answer but squares the condition number, losing roughly twice as many digits.

Time Complexity:
- O(m^2*n) per call to `qr_decompose()`, where $m$ and $n$ are the numbers of rows and columns.
- O(m*n^2) per call to `least_squares()`, since it never forms `q` explicitly.

Space Complexity:
- O(m^2 + m*n) for the matrices returned by `qr_decompose()`.
- O(m*n) auxiliary and O(n) for the vector returned by `least_squares()`.

*/

#include <cassert>
#include <cmath>
#include <optional>
#include <vector>

using Matrix = std::vector<std::vector<double>>;

// Reflects the trailing rows of the given column onto the axis, returning the reflection vector.
std::vector<double> householder_vector(const Matrix &a, int col) {
  int m = static_cast<int>(a.size());
  std::vector<double> v(m);
  double norm = 0;
  for (int i = col; i < m; i++) {
    v[i] = a[i][col];
    norm += v[i] * v[i];
  }
  norm = std::sqrt(norm);
  // Reflect away from a[col][col] so that the subtraction below never cancels.
  v[col] += a[col][col] >= 0 ? norm : -norm;
  return v;
}

void qr_decompose(const Matrix &a, Matrix *q, Matrix *r) {
  int m = static_cast<int>(a.size()), n = a.empty() ? 0 : static_cast<int>(a[0].size());
  assert(m >= n && q != nullptr && r != nullptr);
  *r = a;
  q->assign(m, std::vector<double>(m));
  for (int i = 0; i < m; i++) {
    (*q)[i][i] = 1;
  }
  for (int col = 0; col < n; col++) {
    std::vector<double> v = householder_vector(*r, col);
    double vv = 0;
    for (int i = col; i < m; i++) {
      vv += v[i] * v[i];
    }
    if (vv == 0) {
      continue;  // The column is already zero below the diagonal.
    }
    for (int j = 0; j < n; j++) {  // Apply the reflection to R from the left.
      double dot = 0;
      for (int i = col; i < m; i++) {
        dot += v[i] * (*r)[i][j];
      }
      for (int i = col; i < m; i++) {
        (*r)[i][j] -= 2 * dot / vv * v[i];
      }
    }
    for (int j = 0; j < m; j++) {  // Accumulate the reflection into Q from the right.
      double dot = 0;
      for (int i = col; i < m; i++) {
        dot += v[i] * (*q)[j][i];
      }
      for (int i = col; i < m; i++) {
        (*q)[j][i] -= 2 * dot / vv * v[i];
      }
    }
  }
}

std::optional<std::vector<double>> least_squares(
    const Matrix &a, const std::vector<double> &b, double eps = 1e-10
) {
  int m = static_cast<int>(a.size()), n = a.empty() ? 0 : static_cast<int>(a[0].size());
  assert(m >= n && static_cast<int>(b.size()) == m);
  Matrix r = a;
  std::vector<double> rhs = b;
  for (int col = 0; col < n; col++) {
    std::vector<double> v = householder_vector(r, col);
    double vv = 0;
    for (int i = col; i < m; i++) {
      vv += v[i] * v[i];
    }
    if (vv == 0) {
      continue;
    }
    for (int j = col; j < n; j++) {
      double dot = 0;
      for (int i = col; i < m; i++) {
        dot += v[i] * r[i][j];
      }
      for (int i = col; i < m; i++) {
        r[i][j] -= 2 * dot / vv * v[i];
      }
    }
    double dot = 0;  // The same reflection applied to the right-hand side.
    for (int i = col; i < m; i++) {
      dot += v[i] * rhs[i];
    }
    for (int i = col; i < m; i++) {
      rhs[i] -= 2 * dot / vv * v[i];
    }
  }
  std::vector<double> x(n);
  for (int i = n - 1; i >= 0; i--) {  // Back substitution over the triangular top block.
    if (std::fabs(r[i][i]) <= eps) {
      return std::nullopt;
    }
    double sum = rhs[i];
    for (int j = i + 1; j < n; j++) {
      sum -= r[i][j] * x[j];
    }
    x[i] = sum / r[i][i];
  }
  return x;
}

/*** Example Usage ***/

#include <cmath>
using namespace std;

bool EQ(double a, double b) {
  return fabs(a - b) < 1e-9;
}

int main() {
  Matrix a{{12, -51, 4}, {6, 167, -68}, {-4, 24, -41}};
  Matrix q, r;
  qr_decompose(a, &q, &r);
  for (int i = 0; i < 3; i++) {
    for (int j = 0; j < 3; j++) {
      double product = 0, orthogonal = 0;
      for (int k = 0; k < 3; k++) {
        product += q[i][k] * r[k][j];
        orthogonal += q[k][i] * q[k][j];
      }
      assert(EQ(product, a[i][j]));                // The factors multiply back to a.
      assert(EQ(orthogonal, i == j ? 1.0 : 0.0));  // The columns of q are orthonormal.
      assert(i <= j || EQ(r[i][j], 0.0));          // r is upper triangular.
    }
  }

  // An exactly solvable square system.
  auto exact = least_squares(Matrix{{2, 1}, {1, 3}}, vector<double>{5, 10});
  assert(exact.has_value() && EQ((*exact)[0], 1.0) && EQ((*exact)[1], 3.0));

  // Fit y = c0 + c1*x to four collinear points, then to points with one outlier.
  Matrix design{{1, 0}, {1, 1}, {1, 2}, {1, 3}};
  auto line = least_squares(design, vector<double>{1, 3, 5, 7});
  assert(line.has_value() && EQ((*line)[0], 1.0) && EQ((*line)[1], 2.0));
  auto noisy = least_squares(design, vector<double>{1, 3, 5, 8});
  assert(noisy.has_value() && EQ((*noisy)[0], 0.8) && EQ((*noisy)[1], 2.3));

  // Dependent columns have no unique least-squares solution.
  assert(!least_squares(Matrix{{1, 2}, {2, 4}, {3, 6}}, vector<double>{1, 2, 3}).has_value());
  return 0;
}