1+ # Inputs
2+ # ------
3+ function copy_input (:: typeof (exponential), A:: AbstractMatrix )
4+ return copy! (similar (A, float (eltype (A))), A)
5+ end
6+
7+ copy_input (:: typeof (exponential), A:: Diagonal ) = copy (A)
8+
9+ function check_input (:: typeof (exponential!), A:: AbstractMatrix , expA:: AbstractMatrix , alg:: AbstractAlgorithm )
10+ m = LinearAlgebra. checksquare (A)
11+ @check_size (expA, (m, m))
12+ return @check_scalar (expA, A)
13+ end
14+
15+ function check_input (:: typeof (exponential!), A:: AbstractMatrix , expA:: AbstractMatrix , :: DiagonalAlgorithm )
16+ m = LinearAlgebra. checksquare (A)
17+ isdiag (A) || throw (DimensionMismatch (" diagonal input matrix expected"
18+ @assert expA isa Diagonal
19+ @check_size (expA, (m, m))
20+ @check_scalar (expA, A)
21+ return nothing
22+ end
23+
24+ # Outputs
25+ # -------
26+ initialize_output (:: typeof (exponential!), A:: AbstractMatrix , :: AbstractAlgorithm ) = A
27+
28+ # Implementation
29+ # --------------
30+ function exponential! (A, expA, alg:: MatrixFunctionViaTaylor )
31+ check_input (exponential!, A, expA, alg)
32+ return exponential_via_taylor! (A, expA)
33+ end
34+
35+
36+ module TaylorExponential
37+
38+ """
39+ exponential_via_taylor!(A::Matrix{T}) where T <: AbstractFloat
40+
41+ An implementation of the Fasi & Higham (2018) Taylor-based scaling
42+ and squaring for arbitrary precision.
43+ """
44+ function exponential_via_taylor! (A:: AbstractMatrix{T} , expA:: AbstractMatrix{T} ) where T
45+ n = size (A, 1 )
46+ ϵ = eps (T)
47+ Apowers = fill (A, 8 ) # Preallocate for powers up to A^8, will be resized if needed
48+ for k = 2 : length (Apowers)
49+ Apowers[k] = Apowers[k- 1 ] * A
50+ end
51+ ρA = opnorm (Apowers[end ], 1 )^ (1 / length (Apowers)) # estimate of the spectral radius using Gelfand's formula
52+
53+ # Find m and s such that (ρA/2^s)^(m+1) / (m+1)! < ϵ
54+ m, s = optimal_taylor_order (ρA, ϵ)
55+
56+ # Scale A down by 2^s
57+ A ./= T (2 )^ s
58+
59+ # Evaluate Taylor via Paterson-Stockmeyer approach
60+ X = evaluate_taylor_ps (As, m)
61+
62+ # Squaring to undo the scaling
63+ for _ in 1 : s
64+ X = X * X
65+ end
66+
67+ return X
68+ end
69+
70+ """
71+ Evaluates Taylor series using Paterson-Stockmeyer logic.
72+ """
73+ function evaluate_taylor_ps (A, m)
74+ T = eltype (A)
75+ k = Int (floor (sqrt (m))) # Chunk size for Paterson-Stockmeyer
76+
77+ # Precompute powers A^2, ..., A^k
78+ powers = Vector {typeof(A)} (undef, k)
79+ powers[1 ] = A
80+ for i in 2 : k
81+ powers[i] = powers[i- 1 ] * A
82+ end
83+
84+ # Horner-like evaluation of the outer polynomial
85+ # e^A ≈ ∑ (A^k)^j * P_j(A)
86+ res = zeros (T, size (A))
87+ num_chunks = Int (ceil ((m+ 1 )/ k))
88+
89+ for j in (num_chunks- 1 ): - 1 : 0
90+ # Evaluate sub-polynomial P_j(A) of degree k-1
91+ poly_chunk = zeros (T, size (A))
92+ for i in 0 : k- 1
93+ idx = j* k + i
94+ if idx > m; continue ; end
95+ coeff = 1 / factorial (big (idx))
96+ if i == 0
97+ poly_chunk += T (coeff) * I
98+ else
99+ poly_chunk += T (coeff) * powers[i]
100+ end
101+ end
102+
103+ if j == num_chunks - 1
104+ res = poly_chunk
105+ else
106+ res = res * powers[k] + poly_chunk
107+ end
108+ end
109+ return res
110+ end
111+
112+ function get_ps_costs (max_m) # number of matrix multiplications for an order m polynomial
113+ costs = Dict {Int, Int} ()
114+ for m in 1 : max_m
115+ # Find k in 1:m that minimizes (k-1) + ceil(m/k) - 1
116+ best_c = m # Default naive cost
117+ for k in 1 : m
118+ c = (k - 1 ) + div (m- 1 , k) # == (k - 1) + ceil(Int, m/k) - 1
119+ if c < best_c
120+ best_c = c
121+ end
122+ end
123+ costs[m] = best_c
124+ end
125+ return costs
126+ end
127+
128+ # To filter only for the "efficient" m (where degree increases for the same cost)
129+ function get_efficient_m (max_m:: Int )
130+ costs = get_ps_costs (max_m)
131+ efficient = Dict {Int, Int} ()
132+ next_cost = costs[1 ]
133+ for m in 1 : max_m
134+ cost = next_cost
135+ next_cost = m < max_m ? costs[m+ 1 ] : cost + 1
136+ if cost < next_cost
137+ efficient[m] = cost
138+ end
139+ end
140+ return efficient
141+ end
142+
143+ PS_COSTS = get_efficient_m (100 )
144+
145+ """
146+ Optimizes m and s to minimize cost ≈ m_mults + s.
147+ """
148+ function optimal_taylor_order (ρA, ϵ)
149+ # In a full Fasi-Higham implementation, this would use a small
150+ # search loop or a cost-model lookup.
151+ # Here is a simplified version for BigFloat:
152+ best_m, best_s = 0 , 0
153+ min_cost = typemax (Int)
154+
155+ # Search over efficient Paterson-Stockmeyer degrees
156+ for (m, c) in PS_COSTS
157+ # Calculate s needed for this m: (normA/2^s)^(m+1) / (m+1)! < ϵ
158+ # log2(normA/2^s) * (m+1) - log2((m+1)!) < log2(ϵ)
159+ # log2(normA) - s - [log2((m+1)!) / (m+1)] < log2(ϵ) / (m+1)
160+ term = (log2 (ρA) - (log2 (factorial (big (m+ 1 ))) / (m+ 1 )) - (log2 (ϵ) / (m+ 1 )))
161+ s = max (0 , ceil (Int, term))
162+
163+ cost = c + s
164+ if cost < min_cost || (cost == min_cost && m < best_m)
165+ min_cost = cost
166+ best_m, best_s = m, s
167+ end
168+ end
169+ return best_m, best_s
170+ end
171+
172+ end
0 commit comments