1 /*
2 * Copyright (C) 2023 Alberto Irurueta Carro (alberto@irurueta.com)
3 *
4 * Licensed under the Apache License, Version 2.0 (the "License");
5 * you may not use this file except in compliance with the License.
6 * You may obtain a copy of the License at
7 *
8 * http://www.apache.org/licenses/LICENSE-2.0
9 *
10 * Unless required by applicable law or agreed to in writing, software
11 * distributed under the License is distributed on an "AS IS" BASIS,
12 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
13 * See the License for the specific language governing permissions and
14 * limitations under the License.
15 */
16 package com.irurueta.numerical;
17
18 import com.irurueta.algebra.AlgebraException;
19 import com.irurueta.algebra.Matrix;
20 import com.irurueta.algebra.Utils;
21
22 /**
23 * Estimates exponential of a square matrix.
24 * This is based on Gene H. Golub and Charles F. Van Loan. Matrix Computations. 3rd ed. 1996. p. 572
25 */
26 public class ExponentialMatrixEstimator {
27
28 /**
29 * Default error tolerance of estimated result element-wise.
30 * When this tolerance is used, algorithm achieves an error close to machine precision, but
31 * might be slightly larger than this value (which is about 1.1e-16).
32 */
33 public static final double TOLERANCE = Math.pow(2.0, -53.0);
34
35 /**
36 * Estimates factorial values.
37 */
38 private final DoubleFactorialEstimator factorialEstimator = new DoubleFactorialEstimator();
39
40 /**
41 * Number of rows of matrix to be estimated. Every time the number of rows change, reused
42 * matrices are re-instantiated for computational efficiency.
43 */
44 private int rows = 0;
45
46 /**
47 * Scaled version of provided input matrix. This instance is reused as long as provided input
48 * matrices keep the same size.
49 */
50 private Matrix as;
51
52 /**
53 * Denominator of Padé approximant. This is reused for efficiency while provided input matrices
54 * keep the same size.
55 */
56 private Matrix d;
57
58 /**
59 * Numerator of Padé approximant. This is reused for efficiency while provided input matrices
60 * keep the same size.
61 */
62 private Matrix n;
63
64 /**
65 * Internal matrix reused for efficiency while provided input matrices keep the same size.
66 */
67 private Matrix x;
68
69 /**
70 * Internal matrix reused for efficiency while provided input matrices keep the same size.
71 */
72 private Matrix tmp;
73
74 /**
75 * Internal matrix reused for efficiency while provided input matrices keep the same size.
76 */
77 private Matrix cX;
78
79 /**
80 * Contains copy of estimated result. This is reused for efficiency while provided input
81 * matrices keep the same size.
82 */
83 private Matrix f;
84
85 /**
86 * Estimates exponential of provided matrix with default error tolerance.
87 * When this tolerance is used, algorithm achieves an error close to machine precision, but
88 * might be slightly larger than this value (which is about 1.1e-6).
89 *
90 * @param a matrix to be used for exponential estimation.
91 * @return estimated exponential matrix.
92 * @throws AlgebraException if there are numerical errors.
93 */
94 public Matrix exponential(final Matrix a) throws AlgebraException {
95 return exponential(a, TOLERANCE);
96 }
97
98 /**
99 * Estimates exponential of provided matrix.
100 * Larger tolerance than default one can be used to reduce computational complexity if less
101 * accuracy is required.
102 *
103 * @param a matrix to be used for exponential estimation.
104 * @param tolerance maximum allowed absolute error tolerance element-wise.
105 * @return estimated exponential matrix.
106 * @throws AlgebraException if there are numerical errors.
107 */
108 public Matrix exponential(final Matrix a, final double tolerance) throws AlgebraException {
109 final var result = new Matrix(a.getRows(), a.getColumns());
110 exponential(a, result, tolerance);
111 return result;
112 }
113
114 /**
115 * Estimates exponential of provided matrix with default error tolerance.
116 * When this tolerance is used, algorithm achieves an error close to machine precision, but
117 * might be slightly larger than this value (which is about 1.1e-6).
118 *
119 * @param a matrix to be used for exponential estimation.
120 * @param result instance where result will be stored.
121 * @throws AlgebraException if there are numerical errors.
122 */
123 public void exponential(final Matrix a, final Matrix result) throws AlgebraException {
124 exponential(a, result, TOLERANCE);
125 }
126
127 /**
128 * Estimates exponential of provided matrix.
129 * Larger tolerance than default one can be used to reduce computational complexity if less
130 * accuracy is required.
131 *
132 * @param a matrix to be used for exponential estimation.
133 * @param result instance where result will be stored.
134 * @param tolerance maximum allowed absolute error tolerance element-wise.
135 * @throws AlgebraException if there are numerical errors.
136 */
137 public void exponential(final Matrix a, final Matrix result, final double tolerance) throws AlgebraException {
138 final var aRows = a.getRows();
139
140 if (aRows != a.getColumns()) {
141 throw new IllegalArgumentException("Matrix must be squared");
142 }
143 if (result.getRows() != aRows || result.getColumns() != aRows) {
144 throw new IllegalArgumentException(
145 "Result and provided matrix must have the same size");
146 }
147 if (tolerance < 0.0) {
148 throw new IllegalArgumentException("Tolerance must be zero or greater");
149 }
150
151 final var normMax = normmax(a);
152 final var j = Math.max(0, 1 + (int) Math.floor(Math.log(normMax) / Math.log(2.0)));
153
154 initialize(aRows);
155
156 // scaled version of A
157 as.copyFrom(a);
158 as.multiplyByScalar(1.0 / Math.pow(2.0, j));
159
160 // find order for required tolerance
161 final var q = findQForTolerance(tolerance, normMax);
162
163 var c = 1.0;
164 for (var k = 1; k <= q; k++) {
165 c = c * (q - k + 1) / ((2 * q - k + 1) * k);
166
167 // X = As * X
168 tmp.copyFrom(x);
169 as.multiply(tmp, x);
170
171 // c * X
172 cX.copyFrom(x);
173 cX.multiplyByScalar(c);
174
175 // N = N + c * X
176 n.add(cX);
177
178 // D = D + (-1)^k * cX
179 if (k % 2 != 0) {
180 // k is odd --> (-1)^k = -1
181 cX.multiplyByScalar(-1.0);
182 }
183 // when k is even --> (-1)^k = 1 and there is no sign multiplication needed
184 d.add(cX);
185 }
186
187 // Solve DF = N for F using Gaussian elimination.
188 Utils.solve(d, n, f);
189
190 // now square j times
191 for (var k = 0; k < j; k++) {
192 // F = F^2 = F * F
193 tmp.copyFrom(f);
194 tmp.multiply(f);
195 f.copyFrom(tmp);
196 }
197
198 result.copyFrom(f);
199 }
200
201 /**
202 * Initializes matrices being reused as long as number of rows is preserved for multiple
203 * provided input matrices for efficiency purposes.
204 *
205 * @param rows number of rows.
206 * @throws AlgebraException if an error occurs during instantiation.
207 */
208 private void initialize(final int rows) throws AlgebraException {
209 if (rows != this.rows) {
210 as = new Matrix(rows, rows);
211 d = Matrix.identity(rows, rows);
212 n = Matrix.identity(rows, rows);
213 x = Matrix.identity(rows, rows);
214
215 tmp = new Matrix(rows, rows);
216 cX = new Matrix(rows, rows);
217
218 f = new Matrix(rows, rows);
219
220 this.rows = rows;
221 } else {
222 as.initialize(0.0);
223 Matrix.identity(d);
224 Matrix.identity(n);
225 Matrix.identity(x);
226
227 tmp.initialize(0.0);
228 cX.initialize(0.0);
229
230 f.initialize(0.0);
231 }
232 }
233
234 /**
235 * Estimates relative error achieved by this algorithm for provided input values.
236 *
237 * @param p Padé approximant order of numerator.
238 * @param q Padé approximant order of denominator.
239 * @return estimated relative error.
240 */
241 private double relativeError(final int p, final int q) {
242 final var pFact = factorialEstimator.factorial(p);
243 final var qFact = factorialEstimator.factorial(q);
244 final var pPlusQFact = factorialEstimator.factorial(p + q);
245 final var pPlusQPlusOneFact = factorialEstimator.factorial(p + q + 1);
246
247 return Math.pow(2.0, 3.0 - (p + q)) * pFact / pPlusQFact * qFact / pPlusQPlusOneFact;
248 }
249
250 /**
251 * Finds required order of Padé approximant for provided maximum allowed relative error to be
252 * achieved.
253 *
254 * @param maxRelativeError maximum allowed relative error.
255 * @return Padé approximant order.
256 */
257 private int findQForRelativeError(final double maxRelativeError) {
258 for (var q = 0; ; q++) {
259 final var relativeError = relativeError(q, q);
260 if (relativeError <= maxRelativeError) {
261 return q;
262 }
263 }
264 }
265
266 /**
267 * Finds required order of Padé approximant for provided absolute error tolerance to be
268 * achieved.
269 *
270 * @param tolerance maximum allowed absolute error tolerance.
271 * @param normA infinite norm (maximum absolute value) of provided matrix.
272 * @return Padé approximant order.
273 */
274 private int findQForTolerance(final double tolerance, final double normA) {
275 return findQForRelativeError(tolerance / normA);
276 }
277
278 /**
279 * Estimates infinite norm of provided matrix.
280 * Infinite norm is equivalent to the maximum absolute value of all matrix elements.
281 *
282 * @param a matrix to compute infinite norm for.
283 * @return estimated infinite norm.
284 */
285 private static double normmax(final Matrix a) {
286 var max = 0.0;
287 var buffer = a.getBuffer();
288 for (var v : buffer) {
289 var value = Math.abs(v);
290 if (value > max) {
291 max = value;
292 }
293 }
294 return max;
295 }
296 }