c++ 用 power method 找矩陣兩個 eigenvalue/eigenvector pair

 

#include <math.h>

#include <stdio.h>

#include <vector>

using namespace std;

using F_matrix = vector<vector<double>>;// matrix store raw vectors

void vectorShow(vector<double> v) {

for (int i = 0; i < v.size(); i ++) { printf("%10.4f", v[i]); }

printf("\n");

}

void matrixShow(F_matrix matrix) {

for (int j = 0; j < matrix.size(); j ++) { vectorShow(matrix[j]); }

printf("\n");

}

F_matrix arrayClone(double *array, int w, int h) {// copy array elements one by one

F_matrix temp(h);

for (int i = 0; i < h; i++) {

temp[i].resize(w);

for (int n = 0; n < w; n++) { temp[i][n]= *array ++; }

}

return temp;

}

vector<double> powerEigen(F_matrix &x, int iterations=100, double threshold=1e-16){

int rank = x.size();

if (rank != x[0].size()) return vector<double>(0);

double lambda = 0;

vector<double> y(rank);

vector<double> b(x[0]);// copy 1st vector in x, so that 2nd iteration will have 1st positive value.

for (int k = 0; k < iterations; k ++) {

for (int j = 0; j < rank; j++) { // ŷ = x[j] • b̂

double temp = 0; // Σ x[j][i]*b[i]

for (int i = 0; i < rank; i++) { temp += x[j][i] * b[i]; }

y[j] = temp;

}

double norm2y = 0; // norm2y = ||ŷ||₂

lambda = 0;

for(int i = 0; i < rank; i ++) {

norm2y += y[i] * y[i]; // ||ŷ||₂ = √|ŷᵗŷ| = √Σ(y[i]*y[i])

lambda += b[i] * y[i]; // λ = b•y; => eigenvalue

}

norm2y = sqrt(norm2y);

double dx = 0; // |dx̂|² = Σ(dx[i]²)

for (int i = 0; i < rank; i ++) {// dx̂ = ŷ - λ*b̂;

double temp = y[i] - lambda * b[i];

dx += temp * temp;

}

for (int i = 0; i < rank ; i ++) {// b̂ = ŷ/sqrt(ŷᵗŷ) => eigenvector

b[i] = y[i] / norm2y;

}

if (dx < threshold) break; // dx = sqrt(dx); => no need

}

vector<double> v;// return eigenvalue and eigenvector in one vector

v.push_back(lambda);

for(int i = 0; i < rank ; i ++) { v.push_back(b[i]); }

return v;

}

int main(int argc, char **argv) {

double a[] = {

1, 2, 3, 4,

4, 5, 6, 2,

17, 8, 9, 3,

14, 7, 5, 6

};

F_matrix A = arrayClone(a, 4, 4);

matrixShow(A);

vector<double> v1 = powerEigen(A);// power method to find 1st eigenvalue/eigenvector pair

vectorShow(v1);

double lambda = v1[0];

for (int i = 0; i < A.size(); i++) A[i][i] -= lambda;

vector<double> v2 = powerEigen(A);// power shift method to find 2nd eigenvalue/eigenvector pair

for (int i = 0; i < A.size(); i++) A[i][i] += lambda; // matrix A restore back

v2[0] += v1[0]; // eigenvalue restore back

vectorShow(v2);

}