rust_igraph/algorithms/embedding/
laplacian_spectral_embedding.rs

1//! Laplacian spectral embedding (ALGO-EM-003).
2//!
3//! Computes an `no`-dimensional Euclidean representation of the graph
4//! based on the spectral decomposition of one of several Laplacian
5//! matrices.
6//!
7//! Three Laplacian types are supported for undirected graphs:
8//!
9//! - **D − A** (`LaplacianType::DA`): the combinatorial Laplacian.
10//! - **D⁻¹ᐟ² A D⁻¹ᐟ²** (`LaplacianType::DAD`): the normalized
11//!   adjacency (related to the symmetric normalised Laplacian by
12//!   `L_sym = I − DAD`).
13//! - **I − D⁻¹ᐟ² A D⁻¹ᐟ²** (`LaplacianType::IDAD`): the symmetric
14//!   normalized Laplacian.
15//!
16//! Counterpart of `igraph_laplacian_spectral_embedding()` from
17//! `references/igraph/src/misc/embedding.c:1071`.
18
19use crate::algorithms::community::lanczos::{EigenWhich, lanczos_top_k};
20use crate::algorithms::embedding::adjacency_spectral_embedding::SpectralWhich;
21use crate::core::{Graph, IgraphError, IgraphResult};
22
23/// Which Laplacian variant to use for spectral embedding.
24#[derive(Debug, Clone, Copy, PartialEq, Eq)]
25pub enum LaplacianType {
26    /// D − A: the combinatorial (unnormalized) Laplacian.
27    DA,
28    /// D⁻¹ᐟ² A D⁻¹ᐟ²: the normalized adjacency matrix.
29    DAD,
30    /// I − D⁻¹ᐟ² A D⁻¹ᐟ²: the symmetric normalized Laplacian.
31    IDAD,
32}
33
34/// Result of Laplacian spectral embedding.
35#[derive(Debug, Clone)]
36pub struct LaplacianSpectralEmbeddingResult {
37    /// Eigenvalues in the order selected by `which`.
38    pub eigenvalues: Vec<f64>,
39    /// Embedding matrix: `embedding[v]` is the `no`-dimensional embedding
40    /// of vertex `v`. Each inner vector has length `no`.
41    pub embedding: Vec<Vec<f64>>,
42}
43
44/// Compute the Laplacian spectral embedding of an undirected graph.
45///
46/// Decomposes a Laplacian matrix of the graph into its top eigenpairs
47/// and returns the eigenvectors as vertex embeddings.
48///
49/// `no`: embedding dimension (number of eigenpairs to compute).
50///
51/// `weights`: optional edge weights. Must have length `ecount` if
52/// provided. All weights must be non-negative.
53///
54/// `which`: which eigenvalues to use. See [`SpectralWhich`].
55///
56/// `lap_type`: which Laplacian variant. See [`LaplacianType`].
57///
58/// `scaled`: if true, multiply each eigenvector column by
59/// `sqrt(|eigenvalue|)`, giving `X = U |D|^{1/2}`.
60///
61/// Edge directions are ignored (the graph is treated as undirected).
62///
63/// # Errors
64///
65/// Returns an error if `no` is 0 or exceeds `vcount`, if `weights`
66/// has the wrong length or contains non-finite values, or if any
67/// weight is negative (for `DAD`/`IDAD` types, which require
68/// non-negative weights for the square-root degree normalization).
69///
70/// # Examples
71///
72/// ```
73/// use rust_igraph::{Graph, laplacian_spectral_embedding, SpectralWhich, LaplacianType};
74///
75/// let mut g = Graph::with_vertices(4);
76/// g.add_edge(0, 1).unwrap();
77/// g.add_edge(1, 2).unwrap();
78/// g.add_edge(2, 3).unwrap();
79/// g.add_edge(0, 3).unwrap();
80///
81/// let r = laplacian_spectral_embedding(&g, 2, None, SpectralWhich::LargestAlgebraic,
82///                                       LaplacianType::DA, false).unwrap();
83/// assert_eq!(r.eigenvalues.len(), 2);
84/// assert_eq!(r.embedding.len(), 4);
85/// assert_eq!(r.embedding[0].len(), 2);
86/// ```
87pub fn laplacian_spectral_embedding(
88    graph: &Graph,
89    no: usize,
90    weights: Option<&[f64]>,
91    which: SpectralWhich,
92    lap_type: LaplacianType,
93    scaled: bool,
94) -> IgraphResult<LaplacianSpectralEmbeddingResult> {
95    let n = graph.vcount() as usize;
96
97    if no == 0 {
98        return Err(IgraphError::InvalidArgument(
99            "no must be at least 1".to_string(),
100        ));
101    }
102    if no > n {
103        return Err(IgraphError::InvalidArgument(format!(
104            "no ({no}) exceeds vertex count ({n})"
105        )));
106    }
107
108    if let Some(w) = weights {
109        if w.len() != graph.ecount() {
110            return Err(IgraphError::InvalidArgument(format!(
111                "weights length ({}) differs from edge count ({})",
112                w.len(),
113                graph.ecount()
114            )));
115        }
116        for &wv in w {
117            if !wv.is_finite() {
118                return Err(IgraphError::InvalidArgument(
119                    "edge weights must be finite".to_string(),
120                ));
121            }
122            if wv < 0.0 {
123                return Err(IgraphError::InvalidArgument(
124                    "edge weights must be non-negative for Laplacian embedding".to_string(),
125                ));
126            }
127        }
128    }
129
130    if n == 0 {
131        return Ok(LaplacianSpectralEmbeddingResult {
132            eigenvalues: Vec::new(),
133            embedding: Vec::new(),
134        });
135    }
136
137    if graph.ecount() == 0 {
138        return Ok(LaplacianSpectralEmbeddingResult {
139            eigenvalues: vec![0.0; no],
140            embedding: vec![vec![0.0; no]; n],
141        });
142    }
143
144    let adj = build_adjacency(graph)?;
145    let deg = compute_degree(&adj, weights, n);
146
147    let eigen_which = match which {
148        SpectralWhich::LargestAlgebraic => EigenWhich::LargestAlgebraic,
149        SpectralWhich::SmallestAlgebraic => EigenWhich::SmallestAlgebraic,
150        SpectralWhich::LargestMagnitude => EigenWhich::LargestMagnitude,
151    };
152    let max_iter = n.saturating_mul(50).max(1000);
153
154    let result = match lap_type {
155        LaplacianType::DA => {
156            let mv = |x: &[f64], y: &mut [f64]| matvec_da(&adj, weights, &deg, n, x, y);
157            lanczos_top_k(n, &mv, no, eigen_which, max_iter)
158        }
159        LaplacianType::DAD => {
160            let isd = inv_sqrt_degree(&deg);
161            let mv = |x: &[f64], y: &mut [f64]| matvec_dad(&adj, weights, &isd, n, x, y);
162            lanczos_top_k(n, &mv, no, eigen_which, max_iter)
163        }
164        LaplacianType::IDAD => {
165            let isd = inv_sqrt_degree(&deg);
166            let mv = |x: &[f64], y: &mut [f64]| matvec_idad(&adj, weights, &isd, n, x, y);
167            lanczos_top_k(n, &mv, no, eigen_which, max_iter)
168        }
169    };
170
171    Ok(build_result(result, n, scaled))
172}
173
174// ─── Internal helpers ────────────────────────────────────────────
175
176use crate::algorithms::community::lanczos::LanczosTopKResult;
177
178fn build_result(
179    result: LanczosTopKResult,
180    n: usize,
181    scaled: bool,
182) -> LaplacianSpectralEmbeddingResult {
183    let actual_no = result.eigenvalues.len();
184
185    let mut eigenvalues = result.eigenvalues;
186    zapsmall_vec(&mut eigenvalues);
187
188    let mut embedding = vec![vec![0.0; actual_no]; n];
189    for (d, evec) in result.eigenvectors.iter().enumerate().take(actual_no) {
190        for (v, row) in embedding.iter_mut().enumerate() {
191            row[d] = evec[v];
192        }
193    }
194
195    for row in &mut embedding {
196        zapsmall_vec(row);
197    }
198
199    if scaled {
200        for (d, &eval) in eigenvalues.iter().enumerate().take(actual_no) {
201            let scale = eval.abs().sqrt();
202            for row in &mut embedding {
203                row[d] *= scale;
204            }
205        }
206    }
207
208    LaplacianSpectralEmbeddingResult {
209        eigenvalues,
210        embedding,
211    }
212}
213
214type AdjList = Vec<Vec<(usize, usize)>>; // adj[v] = [(neighbor, edge_id), ...]
215
216fn build_adjacency(graph: &Graph) -> IgraphResult<AdjList> {
217    let n = graph.vcount() as usize;
218    let mut adj: AdjList = vec![Vec::new(); n];
219    for eid in 0..graph.ecount() {
220        #[allow(clippy::cast_possible_truncation)]
221        let eid32 = eid as u32;
222        let s = graph.edge_source(eid32)? as usize;
223        let t = graph.edge_target(eid32)? as usize;
224        adj[s].push((t, eid));
225        if s != t {
226            adj[t].push((s, eid));
227        }
228    }
229    Ok(adj)
230}
231
232fn compute_degree(adj: &AdjList, weights: Option<&[f64]>, n: usize) -> Vec<f64> {
233    let mut deg = vec![0.0; n];
234    for (i, neighbors) in adj.iter().enumerate().take(n) {
235        for &(_nei, eid) in neighbors {
236            let w = match weights {
237                Some(wt) => wt[eid],
238                None => 1.0,
239            };
240            deg[i] += w;
241        }
242    }
243    deg
244}
245
246fn inv_sqrt_degree(deg: &[f64]) -> Vec<f64> {
247    deg.iter()
248        .map(|&d| if d > 0.0 { 1.0 / d.sqrt() } else { 0.0 })
249        .collect()
250}
251
252/// D − A: y = D·x − A·x
253fn matvec_da(
254    adj: &AdjList,
255    weights: Option<&[f64]>,
256    deg: &[f64],
257    n: usize,
258    x: &[f64],
259    y: &mut [f64],
260) {
261    for i in 0..n {
262        y[i] = deg[i] * x[i];
263        for &(nei, eid) in &adj[i] {
264            let w = match weights {
265                Some(wt) => wt[eid],
266                None => 1.0,
267            };
268            y[i] -= w * x[nei];
269        }
270    }
271}
272
273/// D⁻¹ᐟ² A D⁻¹ᐟ²: y = D^{-1/2} (A (D^{-1/2} x))
274fn matvec_dad(
275    adj: &AdjList,
276    weights: Option<&[f64]>,
277    inv_sqrt_deg: &[f64],
278    n: usize,
279    x: &[f64],
280    y: &mut [f64],
281) {
282    // Step 1: z = D^{-1/2} x
283    // Step 2: t = A z
284    // Step 3: y = D^{-1/2} t
285    for i in 0..n {
286        let mut t_i = 0.0;
287        for &(nei, eid) in &adj[i] {
288            let w = match weights {
289                Some(wt) => wt[eid],
290                None => 1.0,
291            };
292            t_i += w * inv_sqrt_deg[nei] * x[nei];
293        }
294        y[i] = inv_sqrt_deg[i] * t_i;
295    }
296}
297
298/// I − D⁻¹ᐟ² A D⁻¹ᐟ²: y = x − D^{-1/2} A D^{-1/2} x
299fn matvec_idad(
300    adj: &AdjList,
301    weights: Option<&[f64]>,
302    inv_sqrt_deg: &[f64],
303    n: usize,
304    x: &[f64],
305    y: &mut [f64],
306) {
307    matvec_dad(adj, weights, inv_sqrt_deg, n, x, y);
308    for i in 0..n {
309        y[i] = x[i] - y[i];
310    }
311}
312
313fn zapsmall_vec(v: &mut [f64]) {
314    if v.is_empty() {
315        return;
316    }
317    let max_abs = v.iter().map(|x| x.abs()).fold(0.0_f64, f64::max);
318    if max_abs == 0.0 {
319        return;
320    }
321    let threshold = max_abs * f64::EPSILON.sqrt();
322    for x in v {
323        if x.abs() < threshold {
324            *x = 0.0;
325        }
326    }
327}
328
329#[cfg(test)]
330mod tests {
331    use super::*;
332
333    fn make_cycle(n: u32) -> Graph {
334        let mut g = Graph::with_vertices(n);
335        for i in 0..n {
336            g.add_edge(i, (i + 1) % n).unwrap();
337        }
338        g
339    }
340
341    fn make_complete(n: u32) -> Graph {
342        let mut g = Graph::with_vertices(n);
343        for i in 0..n {
344            for j in (i + 1)..n {
345                g.add_edge(i, j).unwrap();
346            }
347        }
348        g
349    }
350
351    #[test]
352    fn empty_graph_err() {
353        let g = Graph::with_vertices(0);
354        let r = laplacian_spectral_embedding(
355            &g,
356            1,
357            None,
358            SpectralWhich::LargestAlgebraic,
359            LaplacianType::DA,
360            false,
361        );
362        assert!(r.is_err());
363    }
364
365    #[test]
366    fn no_edges_all_zero() {
367        let g = Graph::with_vertices(4);
368        let r = laplacian_spectral_embedding(
369            &g,
370            2,
371            None,
372            SpectralWhich::LargestAlgebraic,
373            LaplacianType::DA,
374            false,
375        )
376        .unwrap();
377        assert_eq!(r.eigenvalues.len(), 2);
378        for &e in &r.eigenvalues {
379            assert!(e.abs() < 1e-10);
380        }
381    }
382
383    #[test]
384    fn da_complete_k4() {
385        // K4: L = D - A. D = 3I. Eigenvalues of L: 0, 4, 4, 4.
386        let g = make_complete(4);
387        let r = laplacian_spectral_embedding(
388            &g,
389            1,
390            None,
391            SpectralWhich::LargestAlgebraic,
392            LaplacianType::DA,
393            false,
394        )
395        .unwrap();
396        assert!(
397            (r.eigenvalues[0] - 4.0).abs() < 1e-4,
398            "expected 4.0, got {}",
399            r.eigenvalues[0]
400        );
401    }
402
403    #[test]
404    fn da_smallest_is_zero() {
405        // The smallest eigenvalue of the combinatorial Laplacian of a
406        // connected graph is always 0.
407        let g = make_complete(4);
408        let r = laplacian_spectral_embedding(
409            &g,
410            1,
411            None,
412            SpectralWhich::SmallestAlgebraic,
413            LaplacianType::DA,
414            false,
415        )
416        .unwrap();
417        assert!(
418            r.eigenvalues[0].abs() < 1e-4,
419            "expected ~0, got {}",
420            r.eigenvalues[0]
421        );
422    }
423
424    #[test]
425    fn da_path_p4() {
426        // Path P4: L eigenvalues are 0, 2-√2, 2, 2+√2.
427        // Largest = 2 + √2 ≈ 3.414
428        let mut g = Graph::with_vertices(4);
429        g.add_edge(0, 1).unwrap();
430        g.add_edge(1, 2).unwrap();
431        g.add_edge(2, 3).unwrap();
432        let r = laplacian_spectral_embedding(
433            &g,
434            1,
435            None,
436            SpectralWhich::LargestAlgebraic,
437            LaplacianType::DA,
438            false,
439        )
440        .unwrap();
441        let expected = 2.0 + 2.0_f64.sqrt();
442        assert!(
443            (r.eigenvalues[0] - expected).abs() < 1e-3,
444            "expected {expected}, got {}",
445            r.eigenvalues[0]
446        );
447    }
448
449    #[test]
450    fn dad_complete_k4() {
451        // K4: D = 3I, A has eigenvalues 3,-1,-1,-1.
452        // D^{-1/2} A D^{-1/2} = (1/3) A, eigenvalues 1, -1/3, -1/3, -1/3.
453        let g = make_complete(4);
454        let r = laplacian_spectral_embedding(
455            &g,
456            1,
457            None,
458            SpectralWhich::LargestAlgebraic,
459            LaplacianType::DAD,
460            false,
461        )
462        .unwrap();
463        assert!(
464            (r.eigenvalues[0] - 1.0).abs() < 1e-4,
465            "expected 1.0, got {}",
466            r.eigenvalues[0]
467        );
468    }
469
470    #[test]
471    fn idad_complete_k4() {
472        // I - DAD eigenvalues = 1 - DAD eigenvalues = 0, 4/3, 4/3, 4/3.
473        // Largest = 4/3.
474        let g = make_complete(4);
475        let r = laplacian_spectral_embedding(
476            &g,
477            1,
478            None,
479            SpectralWhich::LargestAlgebraic,
480            LaplacianType::IDAD,
481            false,
482        )
483        .unwrap();
484        let expected = 4.0 / 3.0;
485        assert!(
486            (r.eigenvalues[0] - expected).abs() < 1e-4,
487            "expected {expected}, got {}",
488            r.eigenvalues[0]
489        );
490    }
491
492    #[test]
493    fn idad_smallest_is_zero() {
494        // Smallest eigenvalue of L_sym = I - DAD for connected graph is 0.
495        let g = make_cycle(6);
496        let r = laplacian_spectral_embedding(
497            &g,
498            1,
499            None,
500            SpectralWhich::SmallestAlgebraic,
501            LaplacianType::IDAD,
502            false,
503        )
504        .unwrap();
505        assert!(
506            r.eigenvalues[0].abs() < 1e-4,
507            "expected ~0, got {}",
508            r.eigenvalues[0]
509        );
510    }
511
512    #[test]
513    fn da_cycle_c6() {
514        // C6: L eigenvalues are 2 - 2cos(2πk/6) for k=0..5.
515        // = 0, 1, 3, 4, 3, 1. Largest = 4.
516        let g = make_cycle(6);
517        let r = laplacian_spectral_embedding(
518            &g,
519            1,
520            None,
521            SpectralWhich::LargestAlgebraic,
522            LaplacianType::DA,
523            false,
524        )
525        .unwrap();
526        assert!(
527            (r.eigenvalues[0] - 4.0).abs() < 1e-3,
528            "expected 4.0, got {}",
529            r.eigenvalues[0]
530        );
531    }
532
533    #[test]
534    fn scaled_embedding() {
535        let g = make_complete(4);
536        let r_plain = laplacian_spectral_embedding(
537            &g,
538            1,
539            None,
540            SpectralWhich::LargestAlgebraic,
541            LaplacianType::DA,
542            false,
543        )
544        .unwrap();
545        let r_scaled = laplacian_spectral_embedding(
546            &g,
547            1,
548            None,
549            SpectralWhich::LargestAlgebraic,
550            LaplacianType::DA,
551            true,
552        )
553        .unwrap();
554
555        let scale = r_plain.eigenvalues[0].abs().sqrt();
556        for v in 0..4 {
557            let expected = r_plain.embedding[v][0] * scale;
558            assert!(
559                (r_scaled.embedding[v][0] - expected).abs() < 1e-8,
560                "vertex {v}: expected {expected}, got {}",
561                r_scaled.embedding[v][0]
562            );
563        }
564    }
565
566    #[test]
567    fn weighted_da() {
568        // Star S3 with weights 2: center=0, leaves=1,2,3.
569        // D = diag(6,2,2,2), A has edges with weight 2.
570        // L = D - A. Largest eigenvalue of L(star) with weight w:
571        // eigenvalues are 0, w*sqrt(3), ... For unweighted star:
572        // L eigenvalues are 0,1,1,n. For weighted star S3 with w=2:
573        // D=diag(6,2,2,2), eigenvalues of L: 0, 2, 2, 8.
574        let mut g = Graph::with_vertices(4);
575        g.add_edge(0, 1).unwrap();
576        g.add_edge(0, 2).unwrap();
577        g.add_edge(0, 3).unwrap();
578        let weights = vec![2.0, 2.0, 2.0];
579        let r = laplacian_spectral_embedding(
580            &g,
581            1,
582            Some(&weights),
583            SpectralWhich::LargestAlgebraic,
584            LaplacianType::DA,
585            false,
586        )
587        .unwrap();
588        assert!(
589            (r.eigenvalues[0] - 8.0).abs() < 1e-3,
590            "expected 8.0, got {}",
591            r.eigenvalues[0]
592        );
593    }
594
595    #[test]
596    fn embedding_dimension() {
597        let g = make_cycle(6);
598        let r = laplacian_spectral_embedding(
599            &g,
600            3,
601            None,
602            SpectralWhich::LargestAlgebraic,
603            LaplacianType::DA,
604            false,
605        )
606        .unwrap();
607        assert_eq!(r.eigenvalues.len(), 3);
608        assert_eq!(r.embedding.len(), 6);
609        for row in &r.embedding {
610            assert_eq!(row.len(), 3);
611        }
612    }
613
614    #[test]
615    fn invalid_no() {
616        let g = Graph::with_vertices(3);
617        let r = laplacian_spectral_embedding(
618            &g,
619            5,
620            None,
621            SpectralWhich::LargestAlgebraic,
622            LaplacianType::DA,
623            false,
624        );
625        assert!(r.is_err());
626    }
627
628    #[test]
629    fn negative_weight_rejected() {
630        let mut g = Graph::with_vertices(3);
631        g.add_edge(0, 1).unwrap();
632        g.add_edge(1, 2).unwrap();
633        let r = laplacian_spectral_embedding(
634            &g,
635            1,
636            Some(&[1.0, -1.0]),
637            SpectralWhich::LargestAlgebraic,
638            LaplacianType::DA,
639            false,
640        );
641        assert!(r.is_err());
642    }
643
644    #[test]
645    fn dad_regular_graph() {
646        // For a d-regular graph, D = dI, so D^{-1/2} A D^{-1/2} = A/d.
647        // C4 is 2-regular. A eigenvalues: 2, 0, -2, 0.
648        // DAD eigenvalues: 1, 0, -1, 0. Largest = 1.
649        let g = make_cycle(4);
650        let r = laplacian_spectral_embedding(
651            &g,
652            1,
653            None,
654            SpectralWhich::LargestAlgebraic,
655            LaplacianType::DAD,
656            false,
657        )
658        .unwrap();
659        assert!(
660            (r.eigenvalues[0] - 1.0).abs() < 1e-4,
661            "expected 1.0, got {}",
662            r.eigenvalues[0]
663        );
664    }
665}
rust_igraph/algorithms/embedding/laplacian_spectral_embedding.rs

rust_igraph/algorithms/embedding/
laplacian_spectral_embedding.rs
