6.5.2 Row Reduction

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

/*

Converts a matrix to reduced row echelon form using Gaussian elimination to solve a system of linear
equations and compute rank. Each round finds a row with a nonzero entry in the current leading
column, swaps the largest-magnitude candidate into place for partial pivoting, normalizes that row,
and subtracts multiples of it from every other row to clear the column. In practice, this method can
still be prone to rounding error on certain matrices. For a more accurate algorithm for solving
systems of linear equations, LU decomposition with row partial pivoting should be used.

- `row_reduce(a)` assigns the matrix `a` to its reduced row echelon form, returning a reference to
  the modified argument itself.
- `matrix_rank(a)` returns the rank of matrix `a`, i.e. the number of nonzero rows after row
  reduction.
- `solve_system(a, b, &x)` solves the system of linear equations $ax = b$ given an $m$ by $n$ matrix
  `a` of real values, and a length $m$ vector `b`, returning 0 if there is one solution, $-1$
  if there are zero solutions, or $-2$ if there are infinite solutions. If there is exactly one
  solution, then the output pointer `x` is populated with the solution vector of length $n$.

Time Complexity:
- O(m^2*n) per call to `row_reduce(a)`, `matrix_rank(a)`, and `solve_system(a)`, where $m$ and $n$
  are the number of rows and columns of `a`, respectively.

Space Complexity:
- O(1) auxiliary for `row_reduce(a)`.
- O(m*n) auxiliary heap space for `matrix_rank(a)` and `solve_system(a)`.

*/

#include <algorithm>
#include <cmath>
#include <cstddef>
#include <stdexcept>
#include <vector>

const double EPS = 1e-9;

template<class Matrix>
Matrix &row_reduce(Matrix &a) {
  if (a.empty()) {
    return a;
  }
  int rows = static_cast<int>(a.size()), cols = static_cast<int>(a[0].size());
  for (int r = 0, lead = 0; r < rows && lead < cols; lead++) {
    int pivot = r;
    for (int i = r + 1; i < rows; i++) {
      if (fabs(a[i][lead]) > fabs(a[pivot][lead])) {
        pivot = i;
      }
    }
    if (fabs(a[pivot][lead]) < EPS) {
      continue;
    }
    std::swap(a[pivot], a[r]);
    auto lv = a[r][lead];
    for (int j = 0; j < cols; j++) {
      a[r][j] /= lv;
    }
    for (int i = 0; i < rows; i++) {
      if (i != r) {
        lv = a[i][lead];
        for (int j = 0; j < cols; j++) {
          a[i][j] -= lv * a[r][j];
        }
      }
    }
    for (int j = 0; j < lead; j++) {
      a[r][j] = 0;
    }
    a[r][lead] = 1;
    r++;
  }
  return a;
}

template<class Matrix>
int matrix_rank(Matrix a) {
  row_reduce(a);
  int rank = 0;
  for (int i = 0; i < static_cast<int>(a.size()); i++) {
    for (int j = 0; j < static_cast<int>(a[i].size()); j++) {
      if (fabs(a[i][j]) > EPS) {
        rank++;
        break;
      }
    }
  }
  return rank;
}

template<class Matrix, class T>
int solve_system(const Matrix &a, const std::vector<T> &b, std::vector<T> *x) {
  if (x == nullptr || a.empty() || a.size() != b.size()) {
    return -1;
  }
  int rows = static_cast<int>(a.size()), cols = static_cast<int>(a[0].size());
  if (rows < cols) {
    return -2;
  }
  Matrix m(a);
  for (int i = 0; i < rows; i++) {
    m[i].push_back(b[i]);
  }
  row_reduce(m);
  int rank = 0;
  for (int i = 0; i < rows; i++) {
    int lead = -1;
    for (int j = 0; j < cols && lead < 0; j++) {
      if (fabs(m[i][j]) > EPS) {
        lead = j;
      }
    }
    if (lead < 0) {
      // A zero coefficient row with a nonzero constant means 0 = nonzero (no solution).
      if (fabs(m[i][cols]) > EPS) {
        return -1;
      }
    } else {
      rank++;
    }
  }
  // A unique solution requires a pivot in every column; fewer means free variables remain. This
  // also catches trailing free columns, which a per-row `lead > i` test alone would miss.
  if (rank < cols) {
    return -2;
  }
  x->resize(cols);
  for (int i = 0; i < cols; i++) {
    (*x)[i] = m[i][cols];
  }
  return 0;
}

/*** Example Usage ***/

#include <cassert>
using namespace std;

int main() {
  vector<vector<double>> a{{-1, 2, 5}, {1, 0, -6}, {-4, 2, 2}};
  vector<double> b{3, 1, -2};
  vector<double> x;
  assert(solve_system(a, b, &x) == 0);
  assert(matrix_rank(a) == 3);
  assert(matrix_rank(vector<vector<double>>{{1, 2, 3}, {2, 4, 6}, {0, 1, 1}}) == 2);
  for (int i = 0; i < static_cast<int>(a.size()); i++) {
    double sum = 0;
    for (int j = 0; j < static_cast<int>(a[i].size()); j++) {
      sum += a[i][j] * x[j];
    }
    assert(fabs(sum - b[i]) < EPS);
  }
  return 0;
}