Pseudoinverse (Moore-Penrose)

1
#include <bits/stdc++.h>
2
using namespace std;
3
4
// Utility functions for small dense matrices using vector<vector<double>>
5
using Matrix = vector<vector<double>>;
6
7
Matrix zeros(int r, int c) { return Matrix(r, vector<double>(c, 0.0)); }
8
Matrix identity(int n) { Matrix I = zeros(n, n); for (int i=0;i<n;++i) I[i][i]=1.0; return I; }
9
10
Matrix transpose(const Matrix &A){
11
    int m = (int)A.size(), n = (int)A[0].size();
12
    Matrix AT = zeros(n, m);
13
    for(int i=0;i<m;++i) for(int j=0;j<n;++j) AT[j][i] = A[i][j];
14
    return AT;
15
}
16
17
Matrix matmul(const Matrix &A, const Matrix &B){
18
    int m = (int)A.size(), k = (int)A[0].size();
19
    int n = (int)B[0].size();
20
    Matrix C = zeros(m, n);
21
    for(int i=0;i<m;++i){
22
        for(int t=0;t<k;++t){
23
            double a = A[i][t];
24
            if (a==0.0) continue;
25
            for(int j=0;j<n;++j){
26
                C[i][j] += a * B[t][j];
27
            }
28
        }
29
    }
30
    return C;
31
}
32
33
// Cholesky factorization for SPD matrix G: returns lower triangular L such that G = L L^T
34
Matrix cholesky(const Matrix &G){
35
    int n = (int)G.size();
36
    Matrix L = zeros(n, n);
37
    for(int i=0;i<n;++i){
38
        for(int j=0;j<=i;++j){
39
            double sum = G[i][j];
40
            for(int k=0;k<j;++k) sum -= L[i][k]*L[j][k];
41
            if(i==j){
42
                if(sum <= 0) throw runtime_error("Matrix not SPD (Cholesky failed)");
43
                L[i][i] = sqrt(sum);
44
            } else {
45
                L[i][j] = sum / L[j][j];
46
            }
47
        }
48
    }
49
    return L;
50
}
51
52
// Solve L Y = B for Y, where L is lower triangular; B can have multiple RHS (columns)
53
Matrix solveLower(const Matrix &L, const Matrix &B){
54
    int n = (int)L.size();
55
    int m = (int)B[0].size();
56
    Matrix Y = zeros(n, m);
57
    for(int col=0; col<m; ++col){
58
        for(int i=0;i<n;++i){
59
            double sum = B[i][col];
60
            for(int k=0;k<i;++k) sum -= L[i][k]*Y[k][col];
61
            Y[i][col] = sum / L[i][i];
62
        }
63
    }
64
    return Y;
65
}
66
67
// Solve L^T X = Y for X (upper triangular solve)
68
Matrix solveUpperFromLowerT(const Matrix &L, const Matrix &Y){
69
    int n = (int)L.size();
70
    int m = (int)Y[0].size();
71
    Matrix X = zeros(n, m);
72
    for(int col=0; col<m; ++col){
73
        for(int i=n-1;i>=0;--i){
74
            double sum = Y[i][col];
75
            for(int k=i+1;k<n;++k) sum -= L[k][i]*X[k][col]; // using L^T entries
76
            X[i][col] = sum / L[i][i];
77
        }
78
    }
79
    return X;
80
}
81
82
Matrix pseudoinverse_full_column_rank(const Matrix &A){
83
    // A: m x n with full column rank, m >= n
84
    int m = (int)A.size();
85
    int n = (int)A[0].size();
86
    Matrix AT = transpose(A);              // n x m
87
    Matrix G  = matmul(AT, A);             // n x n (Gram matrix)
88
    Matrix L  = cholesky(G);               // n x n lower triangular
89
    // Compute A^+ = (A^T A)^{-1} A^T by solving twice: X = G^{-1} AT = L^{-T} L^{-1} AT
90
    Matrix Y  = solveLower(L, AT);         // solve L Y = A^T
91
    Matrix X  = solveUpperFromLowerT(L, Y);// solve L^T X = Y -> X = A^+
92
    return X; // X is n x m
93
}
94
95
void printMatrix(const Matrix &A, const string &name){
96
    cout << name << " (" << A.size() << "x" << A[0].size() << ")\n";
97
    cout.setf(ios::fixed); cout<<setprecision(6);
98
    for(const auto &row: A){
99
        for(double v: row) cout << setw(12) << v << ' ';
100
        cout << '\n';
101
    }
102
}
103
104
int main(){
105
    // Example: A is 4x2 with full column rank (tall matrix)
106
    Matrix A = {
107
        {1.0, 1.0},
108
        {1.0, 2.0},
109
        {1.0, 3.0},
110
        {1.0, 4.0}
111
    }; // columns: [1,1,1,1]^T and [1,2,3,4]^T
112
113
    Matrix Aplus = pseudoinverse_full_column_rank(A);
114
115
    // Verify properties approximately: A^+ A ≈ I_2 and A A^+ A ≈ A
116
    Matrix ApA = matmul(Aplus, A);       // should be 2x2 identity
117
    Matrix AAp = matmul(A, Aplus);       // projector 4x4
118
    Matrix AApA = matmul(AAp, A);        // should recover A
119
120
    printMatrix(Aplus, "A^+");
121
    printMatrix(ApA,   "A^+ A");
122
    printMatrix(AApA,  "A A^+ A");
123
124
    return 0;
125
}
126

1
#include <bits/stdc++.h>
2
using namespace std;
3
using Matrix = vector<vector<double>>;
4
using Vec = vector<double>;
5
6
Matrix zeros(int r, int c){ return Matrix(r, vector<double>(c, 0.0)); }
7
Matrix transpose(const Matrix &A){ int m=A.size(), n=A[0].size(); Matrix AT=zeros(n,m); for(int i=0;i<m;++i) for(int j=0;j<n;++j) AT[j][i]=A[i][j]; return AT; }
8
9
// Modified Gram-Schmidt QR for m x n with m >= n
10
void qr_mgs(const Matrix &A, Matrix &Q, Matrix &R){
11
    int m = (int)A.size(), n = (int)A[0].size();
12
    Q = zeros(m, n); R = zeros(n, n);
13
    vector<vector<double>> V = A; // copy columns to orthogonalize
14
    for(int j=0;j<n;++j){
15
        // Compute norm of current vector
16
        double norm = 0.0;
17
        for(int i=0;i<m;++i) norm += V[i][j]*V[i][j];
18
        R[j][j] = sqrt(max(0.0, norm));
19
        if(R[j][j] == 0.0) throw runtime_error("Rank deficiency detected in QR");
20
        for(int i=0;i<m;++i) Q[i][j] = V[i][j] / R[j][j];
21
        // Orthogonalize remaining columns against Q[:,j]
22
        for(int k=j+1;k<n;++k){
23
            double r = 0.0;
24
            for(int i=0;i<m;++i) r += Q[i][j] * V[i][k];
25
            R[j][k] = r;
26
            for(int i=0;i<m;++i) V[i][k] -= Q[i][j] * r;
27
        }
28
    }
29
}
30
31
// Solve upper triangular R x = y (n x n)
32
Vec back_substitute(const Matrix &R, const Vec &y){
33
    int n = (int)R.size();
34
    Vec x(n, 0.0);
35
    for(int i=n-1;i>=0;--i){
36
        double sum = y[i];
37
        for(int k=i+1;k<n;++k) sum -= R[i][k]*x[k];
38
        x[i] = sum / R[i][i];
39
    }
40
    return x;
41
}
42
43
// Compute x = argmin ||Ax - b||_2 using QR: R x = Q^T b
44
Vec least_squares_qr(const Matrix &A, const Vec &b){
45
    int m = (int)A.size(), n = (int)A[0].size();
46
    if((int)b.size()!=m) throw runtime_error("Dimension mismatch");
47
    Matrix Q, R; qr_mgs(A, Q, R);
48
    // y = Q^T b
49
    Vec y(n, 0.0);
50
    for(int j=0;j<n;++j){
51
        double dot = 0.0;
52
        for(int i=0;i<m;++i) dot += Q[i][j] * b[i];
53
        y[j] = dot;
54
    }
55
    // Solve R x = y
56
    return back_substitute(R, y);
57
}
58
59
int main(){
60
    // Fit a line y = c0 + c1 * x to noisy data via least squares
61
    vector<double> xs = {0,1,2,3,4,5};
62
    vector<double> ys = {1.1, 2.0, 2.9, 4.2, 5.1, 5.9};
63
    int m = (int)xs.size();
64
    // Build A with columns [1, x]
65
    Matrix A(m, vector<double>(2));
66
    for(int i=0;i<m;++i){ A[i][0]=1.0; A[i][1]=xs[i]; }
67
    // b is observations
68
    vector<double> b = ys;
69
70
    // Solve x = [c0, c1]
71
    vector<double> coeff = least_squares_qr(A, b);
72
73
    cout.setf(ios::fixed); cout<<setprecision(6);
74
    cout << "c0 = " << coeff[0] << ", c1 = " << coeff[1] << "\n";
75
76
    // Predict and print residual norm
77
    double rss = 0.0;
78
    for(int i=0;i<m;++i){
79
        double pred = coeff[0] + coeff[1]*xs[i];
80
        double r = pred - ys[i];
81
        rss += r*r;
82
    }
83
    cout << "Residual sum of squares = " << rss << "\n";
84
    return 0;
85
}
86