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
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
|
export module linalgebra:expm;
import std;
import :error;
import :vector;
import :matrix;
import :lu;
// References used throughout this file:
// Higham — "Functions of Matrices: Theory and Computation" (SIAM 2008)
// Al-Mohy & Higham (2009), SIAM J. Matrix Anal. Appl. 31(3)
export namespace linalgebra {
// Matrix exponential via Padé [13/13] approximant + scaling and squaring.
[[nodiscard]] Matrix expm(const Matrix& A);
} // namespace linalgebra
namespace {
// Padé [13/13] coefficients, Higham (2008) Table 10.4.
constexpr double pade_b[14] = {
64764752532480000.0,
32382376266240000.0,
7771770303897600.0,
1187353796428800.0,
129060195264000.0,
10559470521600.0,
670442572800.0,
33522128640.0,
1323241920.0,
40840800.0,
960960.0,
16380.0,
182.0,
1.0,
};
// theta_13: ||A||_1 threshold below which no scaling is needed.
constexpr double theta_13 = 5.371920351148152;
// 1-norm of a matrix.
double one_norm(const linalgebra::Matrix& M) {
const std::size_t n = M.cols();
const std::size_t m = M.rows();
double result = 0.0;
for (std::size_t j = 0; j < n; ++j) {
double col_sum = 0.0;
for (std::size_t i = 0; i < m; ++i) col_sum += std::abs(M(i, j));
result = std::max(result, col_sum);
}
return result;
}
// Evaluate the Padé [13/13] numerator U and denominator V for matrix B.
// Uses the factored evaluation from Higham Algorithm 10.20.
//
// W1 = b[13]*A6 + b[11]*A4 + b[9]*A2
// W2 = b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*I
// Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2
// Z2 = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*I
// W = A6*W1 + W2
// U = B * W
// V = A6*Z1 + Z2
//
// expm(B) ≈ (V - U)^{-1} * (V + U)
std::pair<linalgebra::Matrix, linalgebra::Matrix>
pade13(const linalgebra::Matrix& B) {
const std::size_t n = B.rows();
const linalgebra::Matrix I = linalgebra::Matrix::identity(n);
const linalgebra::Matrix A2 = B * B;
const linalgebra::Matrix A4 = A2 * A2;
const linalgebra::Matrix A6 = A2 * A4;
// W1 = b[13]*A6 + b[11]*A4 + b[9]*A2
linalgebra::Matrix W1(n, n, 0.0);
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < n; ++j)
W1(i, j) = pade_b[13] * A6(i, j) + pade_b[11] * A4(i, j) + pade_b[9] * A2(i, j);
// W2 = b[7]*A6 + b[5]*A4 + b[3]*A2 + b[1]*I
linalgebra::Matrix W2(n, n, 0.0);
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < n; ++j)
W2(i, j) = pade_b[7] * A6(i, j) + pade_b[5] * A4(i, j)
+ pade_b[3] * A2(i, j) + pade_b[1] * I(i, j);
// Z1 = b[12]*A6 + b[10]*A4 + b[8]*A2
linalgebra::Matrix Z1(n, n, 0.0);
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < n; ++j)
Z1(i, j) = pade_b[12] * A6(i, j) + pade_b[10] * A4(i, j) + pade_b[8] * A2(i, j);
// Z2 = b[6]*A6 + b[4]*A4 + b[2]*A2 + b[0]*I
linalgebra::Matrix Z2(n, n, 0.0);
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < n; ++j)
Z2(i, j) = pade_b[6] * A6(i, j) + pade_b[4] * A4(i, j)
+ pade_b[2] * A2(i, j) + pade_b[0] * I(i, j);
// W = A6*W1 + W2
const linalgebra::Matrix W = A6 * W1 + W2;
// U = B * W (numerator)
const linalgebra::Matrix U = B * W;
// V = A6*Z1 + Z2 (denominator)
const linalgebra::Matrix V = A6 * Z1 + Z2;
return {U, V};
}
} // namespace
namespace linalgebra {
Matrix expm(const Matrix& A) {
if (A.rows() != A.cols()) {
std::ostringstream oss;
oss << "expm requires a square matrix, got " << A.rows() << "x" << A.cols();
throw DimensionMismatchError(oss.str());
}
const std::size_t n = A.rows();
if (n == 0) return Matrix::identity(0);
// Determine scaling factor s such that ||A / 2^s||_1 <= theta_13.
const double norm_A = one_norm(A);
int s = 0;
if (norm_A > theta_13) {
s = static_cast<int>(std::ceil(std::log2(norm_A / theta_13)));
if (s < 0) s = 0;
}
// Scale B = A / 2^s.
const double scale = 1.0 / std::ldexp(1.0, s);
Matrix B(n, n, 0.0);
for (std::size_t i = 0; i < n; ++i)
for (std::size_t j = 0; j < n; ++j)
B(i, j) = A(i, j) * scale;
// Compute Padé [13/13] approximant: R = (V - U)^{-1} * (V + U).
auto [U, V] = pade13(B);
// Numerator = V + U, Denominator = V - U.
Matrix Numerator(n, n, 0.0);
Matrix Denominator(n, n, 0.0);
for (std::size_t i = 0; i < n; ++i) {
for (std::size_t j = 0; j < n; ++j) {
Numerator(i, j) = V(i, j) + U(i, j);
Denominator(i, j) = V(i, j) - U(i, j);
}
}
const LUResult lu_denom = lu_factor(Denominator);
Matrix R(n, n, 0.0);
for (std::size_t j = 0; j < n; ++j) {
Vector col(n);
for (std::size_t i = 0; i < n; ++i) col[i] = Numerator(i, j);
const Vector sol = lu_solve(lu_denom, col);
for (std::size_t i = 0; i < n; ++i) R(i, j) = sol[i];
}
for (int i = 0; i < s; ++i) R = R * R;
return R;
}
} // namespace linalgebra
|