sparse_ir_basis/ir_dlr.rs
1//! Discrete Lehmann representations built from an IR basis
2//!
3//! The DLR itself needs no IR basis (see
4//! [`DiscreteLehmannRepresentation::new`]). This module connects the two:
5//! [`IrBasis`] is a [`Basis`] with a kernel (the IR basis
6//! [`FiniteTempBasis`]), [`DlrFromIr`] builds a DLR whose poles come from such
7//! a basis, and [`IrBasis::dlr_transform`] the linear map between the
8//! coefficients of the basis and of any DLR on the same domain.
9
10use crate::basis::FiniteTempBasis;
11use crate::basis_trait::Basis;
12use crate::dlr::{DiscreteLehmannRepresentation, IrDlrTransform};
13use crate::error::Error;
14use crate::kernel::{CentrosymmKernel, KernelProperties};
15use crate::matrix::Mat;
16use crate::traits::{Statistics, StatisticsType};
17
18/// A basis with a kernel: the IR basis [`FiniteTempBasis`]
19///
20/// The kernel supplies the pole weights `w(β, ω)` of a DLR built from the
21/// basis. A DLR is not an `IrBasis`, so no DLR can be built from a DLR.
22pub trait IrBasis<S: StatisticsType>: Basis<S> {
23 /// Kernel type
24 type Kernel: KernelProperties;
25
26 /// The kernel the basis was built from
27 fn kernel(&self) -> &Self::Kernel;
28
29 /// The linear map between the coefficients of this basis and of `dlr`.
30 ///
31 /// # Errors
32 /// * [`Error::KernelStatisticsMismatch`] if the kernel of this basis does not
33 /// support the statistics
34 /// * [`Error::InvalidParameter`] named `basis` if the β of this basis differs from that
35 /// of `dlr`, or its kernel weight vanishes at a pole where the DLR
36 /// weight does not
37 /// * [`Error::OutOfDomain`] named `poles` for a pole of `dlr` outside
38 /// [-ωmax, ωmax] of this basis
39 fn dlr_transform(&self, dlr: &DiscreteLehmannRepresentation<S>) -> Result<IrDlrTransform, Error>
40 where
41 S: 'static,
42 {
43 if S::STATISTICS == Statistics::Fermionic && self.kernel().ypower() == 1 {
44 return Err(Error::KernelStatisticsMismatch);
45 }
46 let (beta, wmax) = (Basis::beta(self), Basis::wmax(self));
47 let rel = |a: f64, b: f64| (a - b).abs() <= 1e-12 * a.abs().max(b.abs());
48 if !rel(beta, dlr.beta()) {
49 return Err(Error::InvalidParameter {
50 name: "basis",
51 value: format!("beta = {beta:?}"),
52 reason: format!("must equal beta = {:?} of the DLR", dlr.beta()),
53 });
54 }
55 if let Some(&pole) = dlr.poles().iter().find(|p| !(-wmax..=wmax).contains(*p)) {
56 return Err(Error::OutOfDomain {
57 name: "poles",
58 value: pole,
59 domain: (-wmax, wmax),
60 });
61 }
62 let mut ratio = Vec::with_capacity(dlr.poles().len());
63 for (&pole, &w) in dlr.poles().iter().zip(dlr.pole_weights()) {
64 let wir = self.kernel().regularizer::<S>(beta, pole);
65 if w == wir {
66 ratio.push(1.0);
67 } else if wir == 0.0 {
68 return Err(Error::InvalidParameter {
69 name: "basis",
70 value: format!("kernel weight 0 at pole {pole:?}"),
71 reason: "must not vanish where the DLR weight does not".to_string(),
72 });
73 } else {
74 ratio.push(w / wir);
75 }
76 }
77 let v_at_poles =
78 crate::sampling::mat_from_matrix(&Basis::evaluate_omega(self, dlr.poles())?)?;
79 let s = Basis::svals(self);
80 let fitmat = Mat::<f64>::from_fn([Basis::size(self), dlr.poles().len()], |idx| {
81 let (l, i) = (idx[0], idx[1]);
82 -s[l] * v_at_poles[[i, l]] * ratio[i]
83 });
84 Ok(IrDlrTransform::from_matrix(fitmat))
85 }
86}
87
88impl<K, S> IrBasis<S> for FiniteTempBasis<K, S>
89where
90 K: KernelProperties + CentrosymmKernel + Clone + 'static,
91 S: StatisticsType + 'static,
92{
93 type Kernel = K;
94
95 fn kernel(&self) -> &K {
96 FiniteTempBasis::kernel(self)
97 }
98}
99
100/// Constructors of a [`DiscreteLehmannRepresentation`] from an IR basis
101///
102/// The returned DLR carries the [`IrDlrTransform`] of its source basis, so
103/// [`DiscreteLehmannRepresentation::from_ir_nd`] and
104/// [`DiscreteLehmannRepresentation::to_ir_nd`] work on it.
105pub trait DlrFromIr<S: StatisticsType>: Sized {
106 /// Create a DLR from an IR basis with custom poles
107 ///
108 /// The tau-domain pole basis is built from the logistic representation, while
109 /// kernel-specific regularizers are preserved for compatible kernels.
110 ///
111 /// # Arguments
112 /// * `basis` - The IR basis to construct DLR from
113 /// * `poles` - Pole positions on the real-frequency axis, in
114 /// [-ωmax, ωmax] of `basis`
115 ///
116 /// # Errors
117 /// * [`Error::KernelStatisticsMismatch`] if the kernel does not support
118 /// the requested statistics (e.g. `RegularizedBoseKernel` with fermionic
119 /// statistics)
120 /// * [`Error::EmptyInput`] if `poles` is empty
121 /// * [`Error::NonFiniteInput`] if a pole is NaN or infinite, and
122 /// [`Error::OutOfDomain`] if a pole is outside [-ωmax, ωmax] of `basis`,
123 /// both named `poles`, for the first such pole
124 /// * [`Error::NotSupported`] for a bosonic pole at 0 if the kernel has a
125 /// `ypower` other than 0 or 1 (the DLR knows the limit at 0 for those
126 /// only)
127 ///
128 /// Duplicate poles are accepted. They make [`DiscreteLehmannRepresentation::from_ir_nd`]
129 /// ill-conditioned: the coefficients of equal poles are not unique,
130 /// although the round trip through [`DiscreteLehmannRepresentation::to_ir_nd`] still recovers the
131 /// IR coefficients.
132 fn from_ir_with_poles<B: IrBasis<S>>(basis: &B, poles: Vec<f64>) -> Result<Self, Error>;
133
134 /// Create a DLR from an IR basis with its default pole locations
135 ///
136 /// Uses the default omega sampling points from the basis.
137 ///
138 /// # Arguments
139 /// * `basis` - The IR basis to construct DLR from
140 ///
141 /// # Errors
142 /// * [`Error::InsufficientDefaultPoles`] if the basis has fewer default
143 /// poles than functions. This can happen with certain kernel types
144 /// (e.g. `RegularizedBoseKernel`) due to numerical precision limitations
145 /// in root finding.
146 /// * [`Error::KernelStatisticsMismatch`] as in [`Self::from_ir_with_poles`]
147 /// * [`Error::NotSupported`] as in [`Self::from_ir_with_poles`]: for a default
148 /// pole at 0 of a bosonic basis whose kernel has a `ypower` other than
149 /// 0 or 1
150 /// * The errors of
151 /// [`Basis::default_omega_sampling_points`]
152 /// (NotSupported for an SVE with too few singular functions)
153 fn from_ir<B: IrBasis<S>>(basis: &B) -> Result<Self, Error>;
154}
155
156impl<S: StatisticsType + 'static> DlrFromIr<S> for DiscreteLehmannRepresentation<S> {
157 fn from_ir_with_poles<B: IrBasis<S>>(basis: &B, poles: Vec<f64>) -> Result<Self, Error> {
158 // RegularizedBoseKernel (ypower == 1) is meaningful only for bosonic
159 // statistics; its regularizer panics for fermionic input. Reject the
160 // combination before computing anything.
161 if S::STATISTICS == Statistics::Fermionic && basis.kernel().ypower() == 1 {
162 return Err(Error::KernelStatisticsMismatch);
163 }
164 // Without poles the fitting matrix has no columns and the DLR has no
165 // functions.
166 if poles.is_empty() {
167 return Err(Error::EmptyInput { name: "poles" });
168 }
169
170 let beta = Basis::beta(basis);
171 let wmax = Basis::wmax(basis);
172 let accuracy = Basis::accuracy(basis);
173 let kernel_ypower = basis.kernel().ypower();
174
175 // Each pole must be finite and in [-ωmax, ωmax], the domain of V_l.
176 // evaluate_omega below checks the domain of the basis again.
177 for (i, &pole) in poles.iter().enumerate() {
178 if !pole.is_finite() {
179 return Err(Error::NonFiniteInput {
180 name: "poles",
181 index: vec![i],
182 value: pole,
183 });
184 }
185 if !(-wmax..=wmax).contains(&pole) {
186 return Err(Error::OutOfDomain {
187 name: "poles",
188 value: pole,
189 domain: (-wmax, wmax),
190 });
191 }
192 }
193 // A bosonic pole at 0 is evaluated through its finite limit, which is
194 // known for ypower 0 (regularizer tanh(βω/2)) and 1 (regularizer ω)
195 // only; see zero_pole_tau_limit.
196 if S::STATISTICS == Statistics::Bosonic
197 && !(0..=1).contains(&kernel_ypower)
198 && poles.contains(&0.0)
199 {
200 return Err(Error::NotSupported {
201 what: format!(
202 "a bosonic DLR pole at 0 for a kernel with ypower = {kernel_ypower}: \
203 its limit is known for ypower 0 and 1 only"
204 ),
205 });
206 }
207
208 // Compute fitting matrix: fitmat = -s · V(poles)
209 // This transforms DLR coefficients to IR coefficients
210 let v_at_poles = crate::sampling::mat_from_matrix(&Basis::evaluate_omega(basis, &poles)?)?; // shape: [n_poles, basis_size]
211 let s = Basis::svals(basis); // Non-normalized singular values (same as C++)
212
213 let basis_size = Basis::size(basis);
214 let n_poles = poles.len();
215
216 // fitmat[l, i] = -s[l] * V_l(pole[i])
217 // C++: fitmat = (-A_array * s_array.replicate(1, A.cols())).matrix()
218 let fitmat = Mat::<f64>::from_fn([basis_size, n_poles], |idx| {
219 let l = idx[0];
220 let i = idx[1];
221 -s[l] * v_at_poles[[i, l]]
222 });
223
224 let ir = IrDlrTransform::from_matrix(fitmat);
225
226 let regularizers: Vec<f64> = poles
227 .iter()
228 .map(|&pole| basis.kernel().regularizer::<S>(beta, pole))
229 .collect();
230
231 Ok(DiscreteLehmannRepresentation::from_parts(
232 beta,
233 wmax,
234 accuracy,
235 poles,
236 regularizers,
237 kernel_ypower,
238 Some(ir),
239 ))
240 }
241
242 fn from_ir<B: IrBasis<S>>(basis: &B) -> Result<Self, Error> {
243 let poles = Basis::default_omega_sampling_points(basis)?;
244 let basis_size = Basis::size(basis);
245 if basis_size > poles.len() {
246 return Err(Error::InsufficientDefaultPoles {
247 basis_size,
248 n_poles: poles.len(),
249 });
250 }
251 Self::from_ir_with_poles(basis, poles)
252 }
253}