numerics
Loading...
Searching...
No Matches
power.cpp
Go to the documentation of this file.
1/// @file eigen/power.cpp
2/// @brief Power iteration, inverse iteration, Rayleigh quotient iteration
3
4#include "linalg/eigen/power.hpp"
5#include <cmath>
6#include <stdexcept>
7
8namespace num {
9
10
11PowerResult power_iteration(const Matrix& A,
12 real tol, idx max_iter, Backend backend) {
13 constexpr real tiny = 1e-300;
14 const idx n = A.rows();
15 if (A.cols() != n)
16 throw std::invalid_argument("power_iteration: matrix must be square");
17
18 Vector v(n, 0.0);
19 v[0] = 1.0;
20
21 real lambda = 0.0;
22 PowerResult result{0.0, v, 0, false};
23
24 for (idx iter = 0; iter < max_iter; ++iter) {
25 result.iterations = iter + 1;
26
27 Vector w(n);
28 matvec(A, v, w, backend);
29
30 real new_lambda = dot(v, w, backend);
31 detail::normalise(w);
32
33 real delta = std::abs(new_lambda - lambda);
34 lambda = new_lambda;
35 v = w;
36
37 if (delta < tol * (std::abs(lambda) + tiny)) {
38 result.converged = true;
39 break;
40 }
41 }
42
43 result.eigenvalue = lambda;
44 result.eigenvector = v;
45 return result;
46}
47
48PowerResult inverse_iteration(const Matrix& A, real sigma,
49 real tol, idx max_iter, Backend backend) {
50 constexpr real tiny = 1e-300;
51 const idx n = A.rows();
52 if (A.cols() != n)
53 throw std::invalid_argument("inverse_iteration: matrix must be square");
54
55 // Factorize (A - sigma*I) once
56 Matrix M = A;
57 for (idx i = 0; i < n; ++i) M(i,i) -= sigma;
58 LUResult f = lu(M);
59
60 Vector v(n, 0.0);
61 v[0] = 1.0;
62
63 real lambda = 0.0;
64 PowerResult result{0.0, v, 0, false};
65
66 for (idx iter = 0; iter < max_iter; ++iter) {
67 result.iterations = iter + 1;
68
69 Vector w(n);
70 lu_solve(f, v, w);
71 detail::normalise(w);
72
73 // Rayleigh quotient as eigenvalue estimate
74 Vector Av(n);
75 matvec(A, w, Av, backend);
76 real new_lambda = dot(w, Av, backend);
77
78 real delta = std::abs(new_lambda - lambda);
79 lambda = new_lambda;
80 v = w;
81
82 if (delta < tol * (std::abs(lambda) + tiny)) {
83 result.converged = true;
84 break;
85 }
86 }
87
88 result.eigenvalue = lambda;
89 result.eigenvector = v;
90 return result;
91}
92
93PowerResult rayleigh_iteration(const Matrix& A, const Vector& x0,
94 real tol, idx max_iter, Backend backend) {
95 const idx n = A.rows();
96 if (A.cols() != n)
97 throw std::invalid_argument("rayleigh_iteration: matrix must be square");
98 if (x0.size() != n)
99 throw std::invalid_argument("rayleigh_iteration: x0 size mismatch");
100
101 Vector v = x0;
102 detail::normalise(v);
103
104 // Initial Rayleigh quotient
105 Vector Av(n);
106 matvec(A, v, Av, backend);
107 real sigma = dot(v, Av, backend);
108
109 PowerResult result{sigma, v, 0, false};
110
111 for (idx iter = 0; iter < max_iter; ++iter) {
112 result.iterations = iter + 1;
113
114 // Factorize (A - sigma*I); fresh each iteration (cubic convergence)
115 Matrix M = A;
116 for (idx i = 0; i < n; ++i) M(i,i) -= sigma;
117 LUResult f = lu(M);
118
119 if (f.singular) break;
120
121 Vector w(n);
122 lu_solve(f, v, w);
123 detail::normalise(w);
124
125 matvec(A, w, Av, backend);
126 real new_sigma = dot(w, Av, backend);
127
128 // Residual: ||A*v - sigma*v||
129 real res = 0;
130 for (idx i = 0; i < n; ++i) {
131 real d = Av[i] - new_sigma * w[i];
132 res += d*d;
133 }
134 res = std::sqrt(res);
135
136 sigma = new_sigma;
137 v = w;
138
139 if (res < tol) {
140 result.converged = true;
141 break;
142 }
143 }
144
145 result.eigenvalue = sigma;
146 result.eigenvector = v;
147 return result;
148}
149
150} // namespace num
num::Matrix
Dense row-major matrix with optional GPU storage.
Definition matrix.hpp:12
num::detail::normalise
real normalise(Vector &v)
Normalise v in-place; returns the old norm.
Definition power.hpp:33
num::real
double real
Definition types.hpp:10
num::lu
LUResult lu(const Matrix &A)
LU factorization of a square matrix A with partial pivoting.
Definition lu.cpp:21
num::Backend
Backend
Selects which backend handles a linalg operation.
Definition policy.hpp:19
num::ipow
constexpr T ipow(T x) noexcept
Compute x^N at compile time via repeated squaring.
Definition integer_pow.hpp:34
num::idx
std::size_t idx
Definition types.hpp:11
num::matvec
void matvec(const Matrix &A, const Vector &x, Vector &y, Backend b=default_backend)
y = A * x
Definition matrix.cpp:94
num::inverse_iteration
PowerResult inverse_iteration(const Matrix &A, real sigma, real tol=1e-10, idx max_iter=1000, Backend backend=default_backend)
Inverse iteration – finds the eigenvalue closest to a shift sigma.
Definition power.cpp:48
num::dot
real dot(const Vector &x, const Vector &y, Backend b=default_backend)
dot product
Definition vector.cpp:57
num::power_iteration
PowerResult power_iteration(const Matrix &A, real tol=1e-10, idx max_iter=1000, Backend backend=default_backend)
Power iteration – finds the eigenvalue largest in absolute value.
Definition power.cpp:11
num::e
constexpr real e
Definition math.hpp:41
num::lu_solve
void lu_solve(const LUResult &f, const Vector &b, Vector &x)
Solve A*x = b using a precomputed LU factorization.
Definition lu.cpp:67
num::rayleigh_iteration
PowerResult rayleigh_iteration(const Matrix &A, const Vector &x0, real tol=1e-10, idx max_iter=50, Backend backend=default_backend)
Rayleigh quotient iteration – cubically convergent.
Definition power.cpp:93
power.hpp
Power iteration, inverse iteration, Rayleigh quotient iteration.
num::LUResult
Result of an LU factorization with partial pivoting (PA = LU)
Definition lu.hpp:19
num::PowerResult
Result of a single-eigenvalue iteration.
Definition power.hpp:24