tensor4all_core/defaults/
svd.rs

1//! SVD decomposition for tensors.
2//!
3//! Provides [`svd`] and [`svd_with`] for computing truncated SVD of
4//! [`IdxTensor`] values. The tensor is unfolded into a matrix by
5//! splitting its indices into left and right groups, then the standard
6//! matrix SVD is computed and truncated according to [`SvdOptions`].
7//!
8//! This module works with concrete types (`DynIndex`, `IdxTensor`) only.
9
10use crate::defaults::idx_tensor::unfold_split_inner;
11use crate::defaults::DynIndex;
12use crate::index_like::IndexLike;
13use crate::truncation::{
14    validate_svd_truncation_options, validate_svd_truncation_policy, SingularValueMeasure,
15    SvdTruncationOptionsError, SvdTruncationPolicy, ThresholdScale, TruncationRule,
16};
17use crate::{ExecutionContext, IdxTensor};
18use num_complex::Complex64;
19use std::sync::Mutex;
20use tenferro::DType;
21use tenferro_ad::EagerTensor;
22use tenferro_linalg::EagerTensorLinalgExt;
23#[cfg(test)]
24use tensor4all_tensorbackend::native_tensor_primal_to_dense_col_major;
25use thiserror::Error;
26
27/// Error type for SVD operations in tensor4all-linalg.
28#[derive(Debug, Error)]
29pub enum SvdError {
30    /// SVD computation failed.
31    #[error("SVD computation failed: {0}")]
32    ComputationError(#[from] anyhow::Error),
33    /// Invalid truncation threshold value (must be finite and non-negative).
34    #[error("Invalid SVD truncation threshold: {0}. Threshold must be finite and non-negative.")]
35    InvalidThreshold(f64),
36    /// `max_bond_dim` is set to zero.
37    #[error("max_bond_dim must be at least 1 when set")]
38    InvalidMaxBondDim,
39}
40
41/// Options for SVD decomposition with truncation control.
42///
43/// # Examples
44///
45/// ```
46/// use tensor4all_core::svd::{SvdOptions, svd_with};
47/// use tensor4all_core::{DynIndex, SvdTruncationPolicy, IdxTensor};
48///
49/// let i = DynIndex::new_dyn(3);
50/// let j = DynIndex::new_dyn(3);
51/// let data: Vec<f64> = (0..9).map(|x| x as f64).collect();
52/// let tensor = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
53///
54/// let opts = SvdOptions::new().with_policy(SvdTruncationPolicy::new(1e-10));
55/// let (u, s, v) = svd_with::<f64>(&tensor, &[i.clone()], &opts).unwrap();
56///
57/// // U has left index + bond, S is diagonal bond x bond, V has right index + bond
58/// assert_eq!(u.dims()[0], 3);
59/// assert_eq!(s.dims().len(), 2);
60/// ```
61#[derive(Debug, Clone, Copy)]
62pub struct SvdOptions {
63    /// Maximum retained rank after policy-based truncation.
64    pub max_bond_dim: Option<usize>,
65    /// Per-call SVD truncation policy.
66    /// If `None`, the global default policy is used.
67    pub policy: Option<SvdTruncationPolicy>,
68    truncate: bool,
69}
70
71impl Default for SvdOptions {
72    fn default() -> Self {
73        Self::new()
74    }
75}
76
77impl SvdOptions {
78    /// Create new SVD options with no overrides.
79    #[must_use]
80    pub fn new() -> Self {
81        Self {
82            max_bond_dim: None,
83            policy: None,
84            truncate: true,
85        }
86    }
87
88    /// Set the maximum retained rank.
89    #[must_use]
90    pub fn with_max_bond_dim(mut self, max_bond_dim: usize) -> Self {
91        self.max_bond_dim = Some(max_bond_dim);
92        self
93    }
94
95    /// Set the SVD truncation policy override.
96    #[must_use]
97    pub fn with_policy(mut self, policy: SvdTruncationPolicy) -> Self {
98        self.policy = Some(policy);
99        self
100    }
101
102    pub(crate) fn full_rank() -> Self {
103        Self {
104            max_bond_dim: None,
105            policy: None,
106            truncate: false,
107        }
108    }
109}
110
111fn default_policy_guard() -> std::sync::MutexGuard<'static, SvdTruncationPolicy> {
112    match DEFAULT_SVD_TRUNCATION_POLICY.lock() {
113        Ok(guard) => guard,
114        Err(poisoned) => poisoned.into_inner(),
115    }
116}
117
118// Default value: relative per-value threshold 1e-12.
119static DEFAULT_SVD_TRUNCATION_POLICY: Mutex<SvdTruncationPolicy> =
120    Mutex::new(SvdTruncationPolicy::new(1e-12));
121
122/// Get the global default truncation policy for SVD.
123///
124/// The default policy is `SvdTruncationPolicy::new(1e-12)`.
125#[must_use]
126pub fn default_svd_truncation_policy() -> SvdTruncationPolicy {
127    *default_policy_guard()
128}
129
130/// Set the global default truncation policy for SVD.
131///
132/// # Arguments
133/// * `policy` - SVD truncation policy to use when `SvdOptions::policy` is `None`
134///
135/// # Errors
136/// Returns `SvdError::InvalidThreshold` if `policy.threshold` is invalid.
137pub fn set_default_svd_truncation_policy(policy: SvdTruncationPolicy) -> Result<(), SvdError> {
138    validate_svd_truncation_policy(policy).map_err(|e| SvdError::InvalidThreshold(e.0))?;
139    *default_policy_guard() = policy;
140    Ok(())
141}
142
143fn singular_value_measure(value: f64, measure: SingularValueMeasure) -> f64 {
144    match measure {
145        SingularValueMeasure::Value => value,
146        SingularValueMeasure::SquaredValue => value * value,
147    }
148}
149
150/// Compute the retained rank based on an explicit SVD truncation policy.
151pub(crate) fn compute_retained_rank(s_vec: &[f64], policy: &SvdTruncationPolicy) -> usize {
152    if s_vec.is_empty() {
153        return 1;
154    }
155
156    let measured: Vec<f64> = s_vec
157        .iter()
158        .map(|&value| singular_value_measure(value, policy.measure))
159        .collect();
160    if measured.iter().all(|&value| value == 0.0) {
161        return 1;
162    }
163
164    let retained = match (policy.scale, policy.rule) {
165        (ThresholdScale::Relative, TruncationRule::PerValue) => {
166            let reference = measured.iter().copied().fold(0.0_f64, f64::max);
167            measured
168                .iter()
169                .take_while(|&&value| reference > 0.0 && value / reference > policy.threshold)
170                .count()
171        }
172        (ThresholdScale::Absolute, TruncationRule::PerValue) => measured
173            .iter()
174            .take_while(|&&value| value > policy.threshold)
175            .count(),
176        (ThresholdScale::Relative, TruncationRule::DiscardedTailSum) => {
177            let total: f64 = measured.iter().sum();
178            if total == 0.0 {
179                1
180            } else {
181                let mut discarded = 0.0;
182                let mut keep = measured.len();
183                for (i, value) in measured.iter().enumerate().rev() {
184                    if (discarded + value) / total <= policy.threshold {
185                        discarded += value;
186                        keep = i;
187                    } else {
188                        break;
189                    }
190                }
191                keep
192            }
193        }
194        (ThresholdScale::Absolute, TruncationRule::DiscardedTailSum) => {
195            let mut discarded = 0.0;
196            let mut keep = measured.len();
197            for (i, value) in measured.iter().enumerate().rev() {
198                if discarded + value <= policy.threshold {
199                    discarded += value;
200                    keep = i;
201                } else {
202                    break;
203                }
204            }
205            keep
206        }
207    };
208
209    retained.max(1)
210}
211
212#[cfg(test)]
213fn singular_values_from_native(tensor: &tenferro::Tensor) -> Result<Vec<f64>, SvdError> {
214    match tensor.dtype() {
215        DType::F64 => native_tensor_primal_to_dense_col_major::<f64>(tensor)
216            .map_err(|e| SvdError::ComputationError(e.source)),
217        DType::C64 => native_tensor_primal_to_dense_col_major::<Complex64>(tensor)
218            .map(|values| values.into_iter().map(|value| value.re).collect())
219            .map_err(|e| SvdError::ComputationError(e.source)),
220        other => Err(SvdError::ComputationError(anyhow::anyhow!(
221            "native SVD returned unsupported singular-value scalar type {other:?}"
222        ))),
223    }
224}
225
226fn singular_values_from_eager(tensor: &EagerTensor) -> Result<Vec<f64>, SvdError> {
227    let value = tensor
228        .value()
229        .map_err(|source| SvdError::ComputationError(anyhow::Error::new(source)))?;
230    match tensor.dtype() {
231        DType::F64 => value
232            .as_slice::<f64>()
233            .map(|values| values.to_vec())
234            .map_err(|source| SvdError::ComputationError(anyhow::Error::new(source))),
235        DType::C64 => value
236            .as_slice::<Complex64>()
237            .map(|values| values.iter().map(|value| value.re).collect())
238            .map_err(|source| SvdError::ComputationError(anyhow::Error::new(source))),
239        other => Err(SvdError::ComputationError(anyhow::anyhow!(
240            "eager SVD returned unsupported singular-value scalar type {other:?}"
241        ))),
242    }
243}
244
245/// Read the singular-value decision vector of a resident factor.
246///
247/// The factor stays resident; only its values cross the explicit readback
248/// boundary through `context`.
249#[cfg(feature = "tenferro-cuda")]
250fn resident_singular_values(
251    s_inner: &EagerTensor,
252    context: &ExecutionContext,
253) -> Result<Vec<f64>, SvdError> {
254    let indices: Vec<DynIndex> = s_inner
255        .shape()
256        .iter()
257        .map(|&dim| DynIndex::new_dyn(dim))
258        .collect();
259    let decision =
260        IdxTensor::from_inner(indices, s_inner.clone()).map_err(SvdError::ComputationError)?;
261    decision
262        .read_decision_data(context)
263        .map_err(|error| SvdError::ComputationError(anyhow::Error::new(error)))
264}
265
266type SvdTruncatedEagerResult = (
267    EagerTensor,
268    EagerTensor,
269    EagerTensor,
270    Vec<f64>,
271    DynIndex,
272    Vec<DynIndex>,
273    Vec<DynIndex>,
274);
275
276fn svd_truncated_inner(
277    t: &IdxTensor,
278    left_inds: &[DynIndex],
279    options: &SvdOptions,
280) -> Result<SvdTruncatedEagerResult, SvdError> {
281    if t.is_cuda_resident() {
282        return Err(SvdError::ComputationError(anyhow::anyhow!(
283            "CUDA-resident input requires svd_with_in with the owning execution context"
284        )));
285    }
286    svd_truncated_inner_scoped(t, left_inds, options, None)
287}
288
289/// Compute truncated SVD in a caller-owned execution context.
290///
291/// The factorization and truncation slicing execute in `context`. With a CUDA
292/// context only the singular-value decision vector crosses the explicit
293/// readback boundary (via [`IdxTensor::read_decision_data`); with a CPU
294/// context the decision uses the same host read as [`svd_with`].
295///
296/// # Errors
297///
298/// Returns `SvdError` when the tensor does not belong to `context`, when the
299/// indices or options are invalid, or when the factorization or explicit
300/// decision readback fails.
301pub(crate) fn svd_truncated_inner_in(
302    t: &IdxTensor,
303    left_inds: &[DynIndex],
304    options: &SvdOptions,
305    context: &ExecutionContext,
306) -> Result<SvdTruncatedEagerResult, SvdError> {
307    t.validate_context(context)
308        .map_err(|error| SvdError::ComputationError(anyhow::Error::new(error)))?;
309    svd_truncated_inner_scoped(t, left_inds, options, Some(context))
310}
311
312fn svd_truncated_inner_scoped(
313    t: &IdxTensor,
314    left_inds: &[DynIndex],
315    options: &SvdOptions,
316    context: Option<&ExecutionContext>,
317) -> Result<SvdTruncatedEagerResult, SvdError> {
318    let (matrix_inner, _, m, n, left_indices, right_indices) = unfold_split_inner(t, left_inds)
319        .map_err(|e| anyhow::anyhow!("Failed to unfold tensor: {}", e))
320        .map_err(SvdError::ComputationError)?;
321    let k = m.min(n);
322
323    let (mut u_inner, mut s_inner, mut vt_inner) = matrix_inner
324        .svd()
325        .map_err(|e| SvdError::ComputationError(anyhow::anyhow!("{e}")))?;
326    #[cfg(feature = "tenferro-cuda")]
327    let s_full = if let Some(context) = context {
328        if matches!(context, ExecutionContext::Cuda(_)) {
329            resident_singular_values(&s_inner, context)?
330        } else {
331            singular_values_from_eager(&s_inner)?
332        }
333    } else {
334        singular_values_from_eager(&s_inner)?
335    };
336    #[cfg(not(feature = "tenferro-cuda"))]
337    let s_full = {
338        let _ = context;
339        singular_values_from_eager(&s_inner)?
340    };
341    let mut r = if options.truncate {
342        // Shared root guard: reject `max_bond_dim == Some(0)` and invalid
343        // explicit policies before any truncation, so the caller cannot
344        // silently clamp every factorization to rank one.
345        validate_svd_truncation_options(options.max_bond_dim, options.policy).map_err(|error| {
346            match error {
347                SvdTruncationOptionsError::ZeroMaxBondDim => SvdError::InvalidMaxBondDim,
348                SvdTruncationOptionsError::InvalidThreshold(threshold) => {
349                    SvdError::InvalidThreshold(threshold)
350                }
351            }
352        })?;
353        let policy = options.policy.unwrap_or_else(default_svd_truncation_policy);
354        validate_svd_truncation_policy(policy).map_err(|e| SvdError::InvalidThreshold(e.0))?;
355
356        let mut retained = compute_retained_rank(&s_full, &policy);
357        if let Some(max_bond_dim) = options.max_bond_dim {
358            retained = retained.min(max_bond_dim);
359        }
360        retained.max(1)
361    } else {
362        k.max(1)
363    };
364    r = r.min(s_full.len());
365    if r < k {
366        let keep: Vec<usize> = (0..r).collect();
367        u_inner = u_inner
368            .take_axis(1, &keep)
369            .map_err(|e| SvdError::ComputationError(anyhow::anyhow!("{e}")))?;
370        s_inner = s_inner
371            .take_axis(0, &keep)
372            .map_err(|e| SvdError::ComputationError(anyhow::anyhow!("{e}")))?;
373        vt_inner = vt_inner
374            .take_axis(0, &keep)
375            .map_err(|e| SvdError::ComputationError(anyhow::anyhow!("{e}")))?;
376    }
377
378    let bond_index = DynIndex::new_bond(r)
379        .map_err(|e| anyhow::anyhow!("Failed to create Link index: {:?}", e))
380        .map_err(SvdError::ComputationError)?;
381    let singular_values = s_full[..r].to_vec();
382
383    Ok((
384        u_inner,
385        s_inner,
386        vt_inner,
387        singular_values,
388        bond_index,
389        left_indices,
390        right_indices,
391    ))
392}
393
394/// Compute SVD decomposition of a tensor with arbitrary rank, returning (U, S, V).
395///
396/// # Errors
397///
398/// Returns an error when the operation fails (a shape or index mismatch, or
399/// /// a backend failure).
400///
401/// # Examples
402///
403/// ```
404/// use tensor4all_core::{IdxTensor, DynIndex, svd};
405///
406/// // Create a 2x3 matrix (rank-1 outer product: all-ones)
407/// let i = DynIndex::new_dyn(2);
408/// let j = DynIndex::new_dyn(3);
409/// let data = vec![1.0_f64; 6]; // all-ones 2x3 matrix
410/// let t = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
411///
412/// let (u, s, v) = svd::<f64>(&t, &[i.clone()]).unwrap();
413///
414/// // U: shape (left_dim, bond) = (2, bond)
415/// assert_eq!(u.dims()[0], 2);
416/// // V: shape (right_dim, bond) = (3, bond)
417/// assert_eq!(v.dims()[0], 3);
418/// // S is a diagonal matrix (bond × bond)
419/// assert_eq!(s.dims().len(), 2);
420/// ```
421pub fn svd<T>(
422    t: &IdxTensor,
423    left_inds: &[DynIndex],
424) -> Result<(IdxTensor, IdxTensor, IdxTensor), SvdError> {
425    svd_with::<T>(t, left_inds, &SvdOptions::default())
426}
427
428/// Compute SVD decomposition of a tensor with arbitrary rank, returning (U, S, V).
429///
430/// This function allows per-call control of the truncation policy via `SvdOptions`.
431/// If `options.policy` is `None`, it uses the global default policy.
432///
433/// # Errors
434///
435/// Returns an error when the operation fails (a shape or index mismatch, or
436/// /// a backend failure).
437///
438/// # Examples
439///
440/// ```
441/// use tensor4all_core::{DynIndex, IdxTensor};
442/// use tensor4all_core::svd::{SvdOptions, svd_with};
443///
444/// let i = DynIndex::new_dyn(4);
445/// let j = DynIndex::new_dyn(4);
446/// // Rank-1 matrix
447/// let mut data = vec![0.0_f64; 16];
448/// data[0] = 1.0;
449/// let tensor = IdxTensor::from_dense(vec![i.clone(), j.clone()], data).unwrap();
450///
451/// use tensor4all_core::SvdTruncationPolicy;
452///
453/// // Truncate with a relative per-value threshold => rank 1
454/// let opts = SvdOptions::new().with_policy(SvdTruncationPolicy::new(1e-10));
455/// let (u, s, _v) = svd_with::<f64>(&tensor, &[i.clone()], &opts).unwrap();
456/// assert_eq!(s.dims()[0], 1);  // rank-1
457///
458/// // Truncate with max_bond_dim => capped
459/// let opts = SvdOptions::new().with_max_bond_dim(2);
460/// let (_u, s, _v) = svd_with::<f64>(&tensor, &[i.clone()], &opts).unwrap();
461/// assert!(s.dims()[0] <= 2);
462/// ```
463pub fn svd_with<T>(
464    t: &IdxTensor,
465    left_inds: &[DynIndex],
466    options: &SvdOptions,
467) -> Result<(IdxTensor, IdxTensor, IdxTensor), SvdError> {
468    let (u_inner, s_inner, vt_inner, _singular_values, bond_index, left_indices, right_indices) =
469        svd_truncated_inner(t, left_inds, options)?;
470
471    svd_assemble(
472        u_inner,
473        s_inner,
474        vt_inner,
475        bond_index,
476        left_indices,
477        right_indices,
478    )
479}
480
481/// Compute SVD decomposition in a caller-owned execution context.
482///
483/// The factorization, truncation slicing, and result assembly all execute in
484/// `context`. With a CUDA context only the singular-value decision vector
485/// crosses the explicit readback boundary; with a CPU context the decision
486/// uses the same host read as [`svd_with`].
487///
488/// # Arguments
489///
490/// * `t` - Input tensor, which must belong to `context`.
491/// * `left_inds` - Indices to place on the left (row) side of the unfolded matrix.
492/// * `options` - SVD options including truncation policy and bond cap.
493/// * `context` - Caller-owned execution context owning the input and results.
494///
495/// # Examples
496///
497/// ```
498/// use std::sync::Arc;
499/// use tensor4all_core::svd::{SvdOptions, svd_with_in};
500/// use tensor4all_core::{DynIndex, ExecutionContext, IdxTensor, SvdTruncationPolicy};
501/// use tensor4all_tensorbackend::CpuExecutionContext;
502/// use tenferro_cpu::CpuBackend;
503///
504/// let context = ExecutionContext::Cpu(Arc::new(
505///     CpuExecutionContext::from_backend(CpuBackend::new()),
506/// ));
507/// let i = DynIndex::new_dyn(4);
508/// let j = DynIndex::new_dyn(4);
509/// let mut data = vec![0.0_f64; 16];
510/// data[0] = 1.0;
511/// let tensor = IdxTensor::from_dense_in(&context, vec![i.clone(), j.clone()], data)?;
512/// let opts = SvdOptions::new().with_policy(SvdTruncationPolicy::new(1e-10));
513/// let (u, s, _v) = svd_with_in::<f64>(&tensor, &[i.clone()], &opts, &context)?;
514/// assert_eq!(s.dims()[0], 1);
515/// assert_eq!(u.dims()[0], 4);
516/// # Ok::<(), Box<dyn std::error::Error>>(())
517/// ```
518///
519/// # Errors
520///
521/// Returns `SvdError` when the tensor does not belong to `context`, when the
522/// indices or options are invalid, or when the factorization or explicit
523/// decision readback fails.
524pub fn svd_with_in<T>(
525    t: &IdxTensor,
526    left_inds: &[DynIndex],
527    options: &SvdOptions,
528    context: &ExecutionContext,
529) -> Result<(IdxTensor, IdxTensor, IdxTensor), SvdError> {
530    let (u_inner, s_inner, vt_inner, _singular_values, bond_index, left_indices, right_indices) =
531        svd_truncated_inner_in(t, left_inds, options, context)?;
532
533    svd_assemble(
534        u_inner,
535        s_inner,
536        vt_inner,
537        bond_index,
538        left_indices,
539        right_indices,
540    )
541}
542
543/// Wrap the truncated eager factors as [`IdxTensor`]s.
544fn svd_assemble(
545    u_inner: EagerTensor,
546    s_inner: EagerTensor,
547    vt_inner: EagerTensor,
548    bond_index: DynIndex,
549    left_indices: Vec<DynIndex>,
550    right_indices: Vec<DynIndex>,
551) -> Result<(IdxTensor, IdxTensor, IdxTensor), SvdError> {
552    let mut u_indices = left_indices;
553    u_indices.push(bond_index.clone());
554    let u_dims: Vec<usize> = u_indices.iter().map(|idx| idx.dim).collect();
555    let u_reshaped = u_inner.reshape(&u_dims).map_err(|e| {
556        SvdError::ComputationError(anyhow::anyhow!("eager SVD U reshape failed: {e}"))
557    })?;
558    let u = IdxTensor::from_inner(u_indices, u_reshaped).map_err(SvdError::ComputationError)?;
559
560    // S carries a fresh `sim` leg (ITensors convention S: [l, l'], V: l');
561    // the SAME sim instance is shared with V^H so the returned triple
562    // reconstructs under plain contraction:
563    // U [left..., bond] · S [bond, sim] · V^H [sim, right...].
564    let sim_bond_index = bond_index.sim();
565    let s_indices = vec![bond_index.clone(), sim_bond_index.clone()];
566    let s = IdxTensor::from_diag_inner(s_indices, s_inner).map_err(SvdError::ComputationError)?;
567
568    let mut vh_indices = vec![sim_bond_index];
569    vh_indices.extend(right_indices);
570    let vh_dims: Vec<usize> = vh_indices.iter().map(|idx| idx.dim).collect();
571    let vt_reshaped = vt_inner.reshape(&vh_dims).map_err(|e| {
572        SvdError::ComputationError(anyhow::anyhow!("eager SVD V^T reshape failed: {e}"))
573    })?;
574    let vh = IdxTensor::from_inner(vh_indices, vt_reshaped).map_err(SvdError::ComputationError)?;
575    let perm: Vec<usize> = (1..vh.indices.len()).chain(std::iter::once(0)).collect();
576    let v = vh
577        .conj()
578        .permute(&perm)
579        .map_err(|e| SvdError::ComputationError(anyhow::Error::new(e)))?;
580
581    Ok((u, s, v))
582}
583
584/// SVD result for factorization, returning `V^H` directly.
585pub(crate) struct SvdFactorizeResult {
586    pub u: IdxTensor,
587    pub s: IdxTensor,
588    pub vh: IdxTensor,
589    pub bond_index: DynIndex,
590    pub singular_values: Vec<f64>,
591    pub rank: usize,
592}
593
594/// Compute truncated SVD for factorization, returning `V^H` instead of `V`.
595pub(crate) fn svd_for_factorize(
596    t: &IdxTensor,
597    left_inds: &[DynIndex],
598    options: &SvdOptions,
599) -> Result<SvdFactorizeResult, SvdError> {
600    let (u_inner, s_inner, vt_inner, singular_values, bond_index, left_indices, right_indices) =
601        svd_truncated_inner(t, left_inds, options)?;
602    svd_factorize_assemble(
603        u_inner,
604        s_inner,
605        vt_inner,
606        singular_values,
607        bond_index,
608        left_indices,
609        right_indices,
610    )
611}
612
613/// Compute truncated SVD for factorization in a caller-owned context.
614///
615/// Context-scoped counterpart of [`svd_for_factorize`]: the factorization and
616/// truncation slicing execute in `context`, with only the singular-value
617/// decision vector crossing the explicit readback boundary on CUDA.
618///
619/// # Errors
620///
621/// Returns `SvdError` when the tensor does not belong to `context`, when the
622/// indices or options are invalid, or when the factorization or explicit
623/// decision readback fails.
624pub(crate) fn svd_for_factorize_in(
625    t: &IdxTensor,
626    left_inds: &[DynIndex],
627    options: &SvdOptions,
628    context: &ExecutionContext,
629) -> Result<SvdFactorizeResult, SvdError> {
630    let (u_inner, s_inner, vt_inner, singular_values, bond_index, left_indices, right_indices) =
631        svd_truncated_inner_in(t, left_inds, options, context)?;
632    svd_factorize_assemble(
633        u_inner,
634        s_inner,
635        vt_inner,
636        singular_values,
637        bond_index,
638        left_indices,
639        right_indices,
640    )
641}
642
643/// Wrap the truncated eager factors as [`IdxTensor`]s, returning `V^H` directly.
644fn svd_factorize_assemble(
645    u_inner: EagerTensor,
646    s_inner: EagerTensor,
647    vt_inner: EagerTensor,
648    singular_values: Vec<f64>,
649    bond_index: DynIndex,
650    left_indices: Vec<DynIndex>,
651    right_indices: Vec<DynIndex>,
652) -> Result<SvdFactorizeResult, SvdError> {
653    let rank = singular_values.len();
654    let mut u_indices = left_indices;
655    u_indices.push(bond_index.clone());
656    let u_dims: Vec<usize> = u_indices.iter().map(|idx| idx.dim).collect();
657    let u_reshaped = u_inner.reshape(&u_dims).map_err(|e| {
658        SvdError::ComputationError(anyhow::anyhow!("eager SVD U reshape failed: {e}"))
659    })?;
660    let u = IdxTensor::from_inner(u_indices, u_reshaped).map_err(SvdError::ComputationError)?;
661
662    let s_indices = vec![bond_index.clone(), bond_index.sim()];
663    let s = IdxTensor::from_diag_inner(s_indices, s_inner).map_err(SvdError::ComputationError)?;
664
665    let mut vh_indices = vec![bond_index.clone()];
666    vh_indices.extend(right_indices);
667    let vh_dims: Vec<usize> = vh_indices.iter().map(|idx| idx.dim).collect();
668    let vt_reshaped = vt_inner.reshape(&vh_dims).map_err(|e| {
669        SvdError::ComputationError(anyhow::anyhow!("eager SVD V^T reshape failed: {e}"))
670    })?;
671    let vh = IdxTensor::from_inner(vh_indices, vt_reshaped).map_err(SvdError::ComputationError)?;
672
673    Ok(SvdFactorizeResult {
674        u,
675        s,
676        vh,
677        bond_index,
678        singular_values,
679        rank,
680    })
681}
682
683#[cfg(test)]
684mod tests;
tensor4all_core/defaults/svd.rs

tensor4all_core/defaults/
svd.rs
