amari_functional/operator/
compact.rs

1//! Compact and Fredholm operators.
2//!
3//! This module provides traits and types for compact operators and
4//! Fredholm operators, which are central to many results in functional analysis.
5//!
6//! # Mathematical Background
7//!
8//! - **Compact operators** map bounded sets to precompact sets. On infinite-dimensional
9//!   spaces, they are the closure of finite-rank operators.
10//!
11//! - **Fredholm operators** have finite-dimensional kernel and cokernel.
12//!   The Fredholm index is ind(T) = dim(ker T) - dim(coker T).
13
14use crate::error::{FunctionalError, Result};
15use crate::operator::traits::{BoundedOperator, LinearOperator};
16use crate::operator::MatrixOperator;
17use crate::phantom::{Bounded, Compact, Fredholm as FredholmMarker};
18use amari_core::Multivector;
19use std::marker::PhantomData;
20
21/// Trait for compact operators.
22///
23/// A compact operator maps bounded sets to precompact sets.
24/// Equivalently, every bounded sequence has a subsequence whose
25/// image converges.
26pub trait CompactOperator<V, W = V>: BoundedOperator<V, W, Bounded> {
27    /// Check if the operator has finite rank.
28    fn is_finite_rank(&self) -> bool;
29
30    /// Get the rank (dimension of range) if finite.
31    fn rank(&self) -> Option<usize>;
32
33    /// Get the singular values (if computable).
34    ///
35    /// Singular values are eigenvalues of √(T*T).
36    fn singular_values(&self) -> Result<Vec<f64>>;
37}
38
39/// Trait for Fredholm operators.
40///
41/// A Fredholm operator has:
42/// - Finite-dimensional kernel (null space)
43/// - Closed range
44/// - Finite-dimensional cokernel
45pub trait FredholmOperator<V, W = V>: BoundedOperator<V, W, Bounded> {
46    /// Compute the dimension of the kernel.
47    fn kernel_dimension(&self) -> usize;
48
49    /// Compute the dimension of the cokernel.
50    fn cokernel_dimension(&self) -> usize;
51
52    /// Compute the Fredholm index: dim(ker T) - dim(coker T).
53    fn index(&self) -> i64 {
54        self.kernel_dimension() as i64 - self.cokernel_dimension() as i64
55    }
56
57    /// Check if the operator is an isomorphism (index 0, trivial kernel).
58    fn is_isomorphism(&self) -> bool {
59        self.kernel_dimension() == 0 && self.cokernel_dimension() == 0
60    }
61}
62
63/// A finite-rank operator represented as a sum of rank-1 operators.
64///
65/// T = Σᵢ |vᵢ⟩⟨wᵢ|
66///
67/// where |vᵢ⟩⟨wᵢ| maps x to ⟨wᵢ, x⟩vᵢ.
68#[derive(Clone)]
69pub struct FiniteRankOperator<const P: usize, const Q: usize, const R: usize> {
70    /// The "left" vectors vᵢ (range vectors).
71    left_vectors: Vec<Multivector<P, Q, R>>,
72    /// The "right" vectors wᵢ (functionals).
73    right_vectors: Vec<Multivector<P, Q, R>>,
74}
75
76impl<const P: usize, const Q: usize, const R: usize> std::fmt::Debug
77    for FiniteRankOperator<P, Q, R>
78{
79    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
80        f.debug_struct("FiniteRankOperator")
81            .field("rank", &self.left_vectors.len())
82            .field("signature", &(P, Q, R))
83            .finish()
84    }
85}
86
87impl<const P: usize, const Q: usize, const R: usize> FiniteRankOperator<P, Q, R> {
88    /// Create a new finite-rank operator from left and right vectors.
89    pub fn new(
90        left_vectors: Vec<Multivector<P, Q, R>>,
91        right_vectors: Vec<Multivector<P, Q, R>>,
92    ) -> Result<Self> {
93        if left_vectors.len() != right_vectors.len() {
94            return Err(FunctionalError::dimension_mismatch(
95                left_vectors.len(),
96                right_vectors.len(),
97            ));
98        }
99        Ok(Self {
100            left_vectors,
101            right_vectors,
102        })
103    }
104
105    /// Create a rank-1 operator |v⟩⟨w|.
106    pub fn rank_one(v: Multivector<P, Q, R>, w: Multivector<P, Q, R>) -> Self {
107        Self {
108            left_vectors: vec![v],
109            right_vectors: vec![w],
110        }
111    }
112
113    /// Get the rank of the operator.
114    pub fn get_rank(&self) -> usize {
115        self.left_vectors.len()
116    }
117
118    /// Convert to a matrix representation.
119    pub fn to_matrix(&self) -> Result<MatrixOperator<P, Q, R>> {
120        let n = MatrixOperator::<P, Q, R>::DIM;
121        let mut entries = vec![0.0; n * n];
122
123        for (v, w) in self.left_vectors.iter().zip(self.right_vectors.iter()) {
124            let v_coeffs = v.to_vec();
125            let w_coeffs = w.to_vec();
126
127            for i in 0..n {
128                for j in 0..n {
129                    entries[i * n + j] += v_coeffs[i] * w_coeffs[j];
130                }
131            }
132        }
133
134        MatrixOperator::new(entries, n, n)
135    }
136}
137
138impl<const P: usize, const Q: usize, const R: usize> LinearOperator<Multivector<P, Q, R>>
139    for FiniteRankOperator<P, Q, R>
140{
141    fn apply(&self, x: &Multivector<P, Q, R>) -> Result<Multivector<P, Q, R>> {
142        let x_coeffs = x.to_vec();
143        let mut result = Multivector::<P, Q, R>::zero();
144
145        for (v, w) in self.left_vectors.iter().zip(self.right_vectors.iter()) {
146            let w_coeffs = w.to_vec();
147
148            // Compute ⟨w, x⟩
149            let inner_product: f64 = w_coeffs
150                .iter()
151                .zip(x_coeffs.iter())
152                .map(|(a, b)| a * b)
153                .sum();
154
155            // Add ⟨w, x⟩ v
156            result = result.add(&(v * inner_product));
157        }
158
159        Ok(result)
160    }
161
162    fn domain_dimension(&self) -> Option<usize> {
163        Some(MatrixOperator::<P, Q, R>::DIM)
164    }
165
166    fn codomain_dimension(&self) -> Option<usize> {
167        Some(MatrixOperator::<P, Q, R>::DIM)
168    }
169}
170
171impl<const P: usize, const Q: usize, const R: usize>
172    BoundedOperator<Multivector<P, Q, R>, Multivector<P, Q, R>, Bounded>
173    for FiniteRankOperator<P, Q, R>
174{
175    fn operator_norm(&self) -> f64 {
176        // Upper bound: sum of products of norms
177        let mut sum = 0.0;
178        for (v, w) in self.left_vectors.iter().zip(self.right_vectors.iter()) {
179            let v_norm: f64 = v.to_vec().iter().map(|x| x * x).sum::<f64>().sqrt();
180            let w_norm: f64 = w.to_vec().iter().map(|x| x * x).sum::<f64>().sqrt();
181            sum += v_norm * w_norm;
182        }
183        sum
184    }
185}
186
187impl<const P: usize, const Q: usize, const R: usize> CompactOperator<Multivector<P, Q, R>>
188    for FiniteRankOperator<P, Q, R>
189{
190    fn is_finite_rank(&self) -> bool {
191        true
192    }
193
194    fn rank(&self) -> Option<usize> {
195        Some(self.left_vectors.len())
196    }
197
198    fn singular_values(&self) -> Result<Vec<f64>> {
199        // Compute via T*T matrix eigenvalues
200        let matrix = self.to_matrix()?;
201        let t_star_t = matrix.transpose().multiply(&matrix)?;
202
203        let eigenvalues = crate::spectral::compute_eigenvalues(&t_star_t, 100, 1e-10)?;
204
205        // Singular values are square roots of eigenvalues
206        let mut singular_values: Vec<f64> = eigenvalues
207            .iter()
208            .map(|e| e.value.abs().sqrt())
209            .filter(|&s| s > 1e-10)
210            .collect();
211
212        singular_values.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal));
213        Ok(singular_values)
214    }
215}
216
217impl<const P: usize, const Q: usize, const R: usize> FredholmOperator<Multivector<P, Q, R>>
218    for FiniteRankOperator<P, Q, R>
219{
220    fn kernel_dimension(&self) -> usize {
221        // Kernel has dimension n - rank
222        let n = MatrixOperator::<P, Q, R>::DIM;
223        n.saturating_sub(self.left_vectors.len())
224    }
225
226    fn cokernel_dimension(&self) -> usize {
227        // For finite-rank operators, cokernel has dimension n - rank
228        let n = MatrixOperator::<P, Q, R>::DIM;
229        n.saturating_sub(self.left_vectors.len())
230    }
231}
232
233/// Wrapper to mark a matrix operator as compact.
234#[derive(Debug, Clone)]
235pub struct CompactMatrixOperator<const P: usize, const Q: usize, const R: usize> {
236    matrix: MatrixOperator<P, Q, R>,
237    _marker: PhantomData<Compact>,
238}
239
240impl<const P: usize, const Q: usize, const R: usize> CompactMatrixOperator<P, Q, R> {
241    /// Create a compact operator from a matrix.
242    ///
243    /// Note: In finite dimensions, all bounded operators are compact.
244    pub fn new(matrix: MatrixOperator<P, Q, R>) -> Self {
245        Self {
246            matrix,
247            _marker: PhantomData,
248        }
249    }
250
251    /// Get the underlying matrix.
252    pub fn matrix(&self) -> &MatrixOperator<P, Q, R> {
253        &self.matrix
254    }
255}
256
257impl<const P: usize, const Q: usize, const R: usize> LinearOperator<Multivector<P, Q, R>>
258    for CompactMatrixOperator<P, Q, R>
259{
260    fn apply(&self, x: &Multivector<P, Q, R>) -> Result<Multivector<P, Q, R>> {
261        self.matrix.apply(x)
262    }
263
264    fn domain_dimension(&self) -> Option<usize> {
265        self.matrix.domain_dimension()
266    }
267
268    fn codomain_dimension(&self) -> Option<usize> {
269        self.matrix.codomain_dimension()
270    }
271}
272
273impl<const P: usize, const Q: usize, const R: usize>
274    BoundedOperator<Multivector<P, Q, R>, Multivector<P, Q, R>, Bounded>
275    for CompactMatrixOperator<P, Q, R>
276{
277    fn operator_norm(&self) -> f64 {
278        self.matrix.operator_norm()
279    }
280}
281
282impl<const P: usize, const Q: usize, const R: usize> CompactOperator<Multivector<P, Q, R>>
283    for CompactMatrixOperator<P, Q, R>
284{
285    fn is_finite_rank(&self) -> bool {
286        // All finite-dimensional operators have finite rank
287        true
288    }
289
290    fn rank(&self) -> Option<usize> {
291        // Compute rank via singular value decomposition
292        let singular_values = self.singular_values().ok()?;
293        Some(singular_values.len())
294    }
295
296    fn singular_values(&self) -> Result<Vec<f64>> {
297        let t_star_t = self.matrix.transpose().multiply(&self.matrix)?;
298        let eigenvalues = crate::spectral::compute_eigenvalues(&t_star_t, 100, 1e-10)?;
299
300        let mut singular_values: Vec<f64> = eigenvalues
301            .iter()
302            .map(|e| e.value.abs().sqrt())
303            .filter(|&s| s > 1e-10)
304            .collect();
305
306        singular_values.sort_by(|a, b| b.partial_cmp(a).unwrap_or(std::cmp::Ordering::Equal));
307        Ok(singular_values)
308    }
309}
310
311/// Wrapper to mark a matrix operator as Fredholm.
312#[derive(Debug, Clone)]
313pub struct FredholmMatrixOperator<const P: usize, const Q: usize, const R: usize> {
314    matrix: MatrixOperator<P, Q, R>,
315    kernel_dim: usize,
316    cokernel_dim: usize,
317    _marker: PhantomData<FredholmMarker>,
318}
319
320impl<const P: usize, const Q: usize, const R: usize> FredholmMatrixOperator<P, Q, R> {
321    /// Create a Fredholm operator from a matrix.
322    ///
323    /// Computes the kernel and cokernel dimensions.
324    pub fn new(matrix: MatrixOperator<P, Q, R>) -> Result<Self> {
325        // Compute kernel dimension by finding null space
326        let eigenvalues = crate::spectral::compute_eigenvalues(&matrix, 100, 1e-10)?;
327        let kernel_dim = eigenvalues.iter().filter(|e| e.value.abs() < 1e-10).count();
328
329        // For square matrices, cokernel = kernel dimension
330        let cokernel_dim = kernel_dim;
331
332        Ok(Self {
333            matrix,
334            kernel_dim,
335            cokernel_dim,
336            _marker: PhantomData,
337        })
338    }
339
340    /// Get the underlying matrix.
341    pub fn matrix(&self) -> &MatrixOperator<P, Q, R> {
342        &self.matrix
343    }
344}
345
346impl<const P: usize, const Q: usize, const R: usize> LinearOperator<Multivector<P, Q, R>>
347    for FredholmMatrixOperator<P, Q, R>
348{
349    fn apply(&self, x: &Multivector<P, Q, R>) -> Result<Multivector<P, Q, R>> {
350        self.matrix.apply(x)
351    }
352
353    fn domain_dimension(&self) -> Option<usize> {
354        self.matrix.domain_dimension()
355    }
356
357    fn codomain_dimension(&self) -> Option<usize> {
358        self.matrix.codomain_dimension()
359    }
360}
361
362impl<const P: usize, const Q: usize, const R: usize>
363    BoundedOperator<Multivector<P, Q, R>, Multivector<P, Q, R>, Bounded>
364    for FredholmMatrixOperator<P, Q, R>
365{
366    fn operator_norm(&self) -> f64 {
367        self.matrix.operator_norm()
368    }
369}
370
371impl<const P: usize, const Q: usize, const R: usize> FredholmOperator<Multivector<P, Q, R>>
372    for FredholmMatrixOperator<P, Q, R>
373{
374    fn kernel_dimension(&self) -> usize {
375        self.kernel_dim
376    }
377
378    fn cokernel_dimension(&self) -> usize {
379        self.cokernel_dim
380    }
381}
382
383#[cfg(test)]
384mod tests {
385    use super::*;
386
387    #[test]
388    fn test_rank_one_operator() {
389        let v = Multivector::<2, 0, 0>::from_slice(&[1.0, 0.0, 0.0, 0.0]);
390        let w = Multivector::<2, 0, 0>::from_slice(&[0.0, 1.0, 0.0, 0.0]);
391
392        let op = FiniteRankOperator::rank_one(v, w);
393
394        // Apply to x = (0, 1, 0, 0)
395        // Result should be ⟨w, x⟩ v = 1 * (1, 0, 0, 0) = (1, 0, 0, 0)
396        let x = Multivector::<2, 0, 0>::from_slice(&[0.0, 1.0, 0.0, 0.0]);
397        let result = op.apply(&x).unwrap();
398
399        assert!((result.to_vec()[0] - 1.0).abs() < 1e-10);
400        assert!(result.to_vec()[1].abs() < 1e-10);
401    }
402
403    #[test]
404    fn test_finite_rank_is_compact() {
405        let v = Multivector::<2, 0, 0>::from_slice(&[1.0, 0.0, 0.0, 0.0]);
406        let w = Multivector::<2, 0, 0>::from_slice(&[1.0, 0.0, 0.0, 0.0]);
407
408        let op = FiniteRankOperator::rank_one(v, w);
409
410        assert!(op.is_finite_rank());
411        assert_eq!(op.rank(), Some(1));
412    }
413
414    #[test]
415    fn test_finite_rank_to_matrix() {
416        let v = Multivector::<2, 0, 0>::from_slice(&[1.0, 0.0, 0.0, 0.0]);
417        let w = Multivector::<2, 0, 0>::from_slice(&[1.0, 0.0, 0.0, 0.0]);
418
419        let op = FiniteRankOperator::rank_one(v, w);
420        let matrix = op.to_matrix().unwrap();
421
422        // The matrix should be a projection onto the first coordinate
423        assert!((matrix.get(0, 0) - 1.0).abs() < 1e-10);
424        assert!(matrix.get(0, 1).abs() < 1e-10);
425        assert!(matrix.get(1, 0).abs() < 1e-10);
426    }
427
428    #[test]
429    fn test_compact_matrix_operator() {
430        let id: MatrixOperator<2, 0, 0> = MatrixOperator::identity();
431        let compact = CompactMatrixOperator::new(id);
432
433        assert!(compact.is_finite_rank());
434        assert_eq!(compact.rank(), Some(4));
435    }
436
437    #[test]
438    fn test_fredholm_identity() {
439        let id: MatrixOperator<2, 0, 0> = MatrixOperator::identity();
440        let fredholm = FredholmMatrixOperator::new(id).unwrap();
441
442        // Identity has trivial kernel and cokernel
443        assert_eq!(fredholm.kernel_dimension(), 0);
444        assert_eq!(fredholm.cokernel_dimension(), 0);
445        assert_eq!(fredholm.index(), 0);
446        assert!(fredholm.is_isomorphism());
447    }
448
449    #[test]
450    fn test_fredholm_singular() {
451        // Create a singular matrix (projection onto first 2 dimensions)
452        let proj: MatrixOperator<2, 0, 0> =
453            MatrixOperator::diagonal(&[1.0, 1.0, 0.0, 0.0]).unwrap();
454        let fredholm = FredholmMatrixOperator::new(proj).unwrap();
455
456        // Should have 2-dimensional kernel
457        assert_eq!(fredholm.kernel_dimension(), 2);
458        assert_eq!(fredholm.index(), 0);
459        assert!(!fredholm.is_isomorphism());
460    }
461
462    #[test]
463    fn test_singular_values() {
464        let id: MatrixOperator<2, 0, 0> = MatrixOperator::identity();
465        let compact = CompactMatrixOperator::new(id);
466
467        let sv = compact.singular_values().unwrap();
468
469        // All singular values of identity are 1
470        for s in &sv {
471            assert!((s - 1.0).abs() < 0.1);
472        }
473    }
474}
amari_functional/operator/compact.rs

amari_functional/operator/
compact.rs
