Suppose you have a matrix $A \in \mathbb{R}^{n \times n}$ and you want to invert it. I've already explained how to invert a matrix (English explanation), but I didn't provide any code and / or runtime analysis.
C++ Code ¶
#include <iostream>
#include <cmath>
#include <vector>
using namespace std;
void print(vector< vector<double> > A) {
int n = A.size();
for (int i=0; i<n; i++) {
for (int j=0; j<2*n; j++) {
cout << A[i][j] << "\t";
if (j == n-1) {
cout << "| ";
}
}
cout << "\n";
}
cout << endl;
}
void calculateInverse(vector< vector<double> >& A) {
int n = A.size();
for (int i=0; i<n; i++) {
// Search for maximum in this column
double maxEl = abs(A[i][i]);
int maxRow = i;
for (int k=i+1; k<n; k++) {
if (abs(A[k][i]) > maxEl) {
maxEl = A[k][i];
maxRow = k;
}
}
// Swap maximum row with current row (column by column)
for (int k=i; k<2*n;k++) {
double tmp = A[maxRow][k];
A[maxRow][k] = A[i][k];
A[i][k] = tmp;
}
// Make all rows below this one 0 in current column
for (int k=i+1; k<n; k++) {
double c = -A[k][i]/A[i][i];
for (int j=i; j<2*n; j++) {
if (i==j) {
A[k][j] = 0;
} else {
A[k][j] += c * A[i][j];
}
}
}
}
// Solve equation Ax=b for an upper triangular matrix A
for (int i=n-1; i>=0; i--) {
for (int k=n; k<2*n;k++) {
A[i][k] /= A[i][i];
}
// this is not necessary, but the output looks nicer:
A[i][i] = 1;
for (int rowModify=i-1;rowModify>=0; rowModify--) {
for (int columModify=n;columModify<2*n;columModify++) {
A[rowModify][columModify] -= A[i][columModify]
* A[rowModify][i];
}
// this is not necessary, but the output looks nicer:
A[rowModify][i] = 0;
}
}
}
int main() {
int n;
cin >> n;
vector<double> line(2*n,0);
vector< vector<double> > A(n,line);
// Read input data
for (int i=0; i<n; i++) {
for (int j=0; j<n; j++) {
cin >> A[i][j];
}
}
for (int i=0; i<n; i++) {
A[i][n+i] = 1;
}
// Print input
print(A);
// Calculate solution
calculateInverse(A);
// Print result
cout << "Result:" << endl;
print(A);
}
This code is VERY similar to the code of Gaussian elimination. In fact, I've only changed all occurrences of n+1 to 2*n. I also had to change lines 57-70, as we need to do all operations now on a matrix instead of a vector.
Time complexity ¶
Lines 43-52:
$\displaystyle Operations = \sum_{i=0}^{n-1} \left (\sum_{k=i+1}^{n-1} (\sum_{j=i}^{2n} 1) \right ) = \frac{5}{6} n^3 - \frac{5}{6} n$ (see Wolfram|Alpha)
Lines 63-70:
$\displaystyle Operations = \sum_{i=0}^{n-1} \left (\sum_{k=0}^{i-1} (\sum_{j=n}^{2n} 1) \right ) = \frac{1}{2} n^3 - \frac{1}{2} n$ (see Wolfram|Alpha)
So we need about $\frac{4}{3} n^3 + \mathcal{O}(n^2)$ operations to invert a matrix with Gauß-Elimination.
Space complexity ¶
My algorithm needs space for the inverse matrix, so it is in $\mathcal{O}(n^2)$.
Inverting an upper triangular matrix ¶
Suppose you have an upper triangular matrix $A \in \mathbb{R^{n \times n}}$ that you would like to invert. It could look like this:
$$\begin{pmatrix} 2 & 7 & 1 & 8 & 2\ 0 & 8 & 1 & 8 & 2\ 0 & 0 & 8 & 4 & 5\ 0 & 0 & 0 & 9 & 0\ 0 & 0 & 0 & 0 & 4 \end{pmatrix}$$
How could we improve the algorithm from above to get speed it up?
Well, we don't need lines 24-53 any more, as those lines bring $A$ to an upper triangular form. But we still need lines 63-70. So we can improve the algorithm to a complexity of $\frac{1}{2} n^3 + \mathcal{O}(n^2)$.