<p>So there were multiple issues with my use of ARPACK. Some of them:</p>
<p>Linear algebra structs not having continous memory alignment, improper use of the shift invert mode, not having enough lanczos vectors...</p>
<p>This is a result after having fixed them:</p>
<p><a href="https://i.sstatic.net/pzK79tAf.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/pzK79tAf.png" alt="enter image description here" /></a></p>
<p>In more detail:</p>
<p>Faer's matrices are not guaranteed to be stored in contiguous memory, so using them for high level storage of the results from ARPACK required some care to avoid unaligning the memory.</p>
<p>For the shift invert mode, I was not setting the correct code to the iparam routine.</p>
<p>And I had neglected to read the part of the documentation that suggests having about 3x more lanczos basis vectors than requested eigenvalues.</p>

CC BY-SA 4.0
2024-07-26T22:06:59.403Z
It may be useful to others working with rust if you provide more details on how you solved the three problems you mentioned.
lightxbulb
CC BY-SA 4.0
2024-07-27T16:45:28.527Z
@lightxbulb added
Makogan
CC BY-SA 4.0
2024-07-26T00:55:51.640Z
2024-07-27T16:45:20.467Z
2024-07-27T16:45:20.467Z
Makogan
<p>I am trying to implement <a href="https://graphics.stanford.edu/courses/cs233-21-spring/ReferencedPapers/understand_geometry_01631196.pdf" rel="nofollow noreferrer">manifold harmonics</a>.</p>
<p>To see if I am doing things right I am visualizing the <a href="https://en.wikipedia.org/wiki/Algebraic_connectivity" rel="nofollow noreferrer">Fiedler eigenvector</a>, which is the eigenvector associated with the second smallest eigenvalue.</p>
<p>On my own code base (which is using ARPACK to compute the eigendecomposition), I am getting strange results.</p>
<p>Normally, if I am requesting the k smallest eigenvalues, according to the wikipedia page for a graph laplacian, the first eigenvalue should be 0. However, in my case I am getting these eigenvalues:</p>
<pre class="lang-rs prettyprint-override"><code>[-4211.039032568338],
[2364.2174825349684],
[3668.6929511652374],
[11190.921998283911],
[17989.923028566314],
</code></pre>
<p>The visualization of the eigenvectors is also weird for example I am getting this image for the second eigenvector (sorted from smallest eigenvalue to largest).</p>
<p><a href="https://i.sstatic.net/TMbDbleJ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/TMbDbleJ.png" alt="enter image description here" /></a></p>
<p>I am not sure if it is the input matrix or my use of arpack that is breaking, I would be more inclined of thinking it is arpack because i have tested the matrix in the past.</p>
<p>To try to debug I build geogram to visualize their results for the same mesh. But it doesn't seem to do anything useful.</p>
<p>I loaded the mesh:</p>
<p><a href="https://i.sstatic.net/LcXILBdr.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/LcXILBdr.png" alt="enter image description here" /></a></p>
<p>Clicked on manifold harmonics:</p>
<p><a href="https://i.sstatic.net/4r5p8TLj.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/4r5p8TLj.png" alt="enter image description here" /></a></p>
<p>Applied them:</p>
<p><a href="https://i.sstatic.net/w6mI16Y8.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/w6mI16Y8.png" alt="enter image description here" /></a></p>
<p>For every attribute labeled MH I am geting that the range is 0</p>
<p><a href="https://i.sstatic.net/AcVKrx8J.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/AcVKrx8J.png" alt="enter image description here" /></a></p>
<p>As if the eigenvectors are always 0. I have tried multiple meshes and they all return this result. I was hoping someone could help me.</p>
<p>As suggested in the comments, I am attempting to use shift-invert:</p>
<pre class="lang-rs prettyprint-override"><code>pub fn eigen_decompose_sparse_matrix(
    mat: &amp;SparseColMat&lt;usize, f64&gt;,
) -&gt; (Col&lt;f64&gt;, Mat&lt;f64&gt;)
{
    debug_assert!(mat.nrows() == mat.ncols());
    faer_arpack_symmetric_f64(
        &amp;mat,
        mat.nrows(),
        Eigenvaluekind::LargestMagnitude.as_str(),
        5,
        10,
        3000,
        true,
    )
}

pub enum Eigenvaluekind
{
    LargestMagnitude,
    SmallestMagnitude,
    LargestRealPart,
    SmallestRealPart,
    LargestImaginaryPart,
    SmallestImaginaryPart,
}

impl Eigenvaluekind
{
    fn as_str(&amp;self) -&gt; &amp;'static str
    {
        match self
        {
            Eigenvaluekind::LargestMagnitude =&gt; &quot;LM&quot;,
            Eigenvaluekind::SmallestMagnitude =&gt; &quot;SM&quot;,
            Eigenvaluekind::LargestRealPart =&gt; &quot;LR&quot;,
            Eigenvaluekind::SmallestRealPart =&gt; &quot;SR&quot;,
            Eigenvaluekind::LargestImaginaryPart =&gt; &quot;LI&quot;,
            Eigenvaluekind::SmallestImaginaryPart =&gt; &quot;SI&quot;,
        }
    }
}

// http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf
fn faer_arpack_symmetric_f64(
    mut mat: &amp;SparseColMat&lt;usize, f64&gt;,
    dimension: usize,
    eigenvalue_kind: &amp;str,
    number_of_eigenvalues: usize,
    number_of_basis_vectors: usize,
    maxiter: usize,
    vectors: bool,
) -&gt; (Col&lt;f64&gt;, Mat&lt;f64&gt;)
{
    let mut ido = 0;
    let mut residual: Col&lt;f64&gt; = Col::zeros(dimension);
    let mut arnoldi_vectors = Mat::&lt;f64&gt;::zeros(dimension, number_of_basis_vectors);
    let mut iparam = [0; 11];
    iparam[0] = 1;
    iparam[2] = maxiter as i32;
    iparam[6] = 1;
    let mut ipntr = [0; 14];
    let mut workd = Col::zeros(3 * dimension);
    let lworkl = number_of_basis_vectors * (number_of_basis_vectors + 8);
    let mut workl = Col::&lt;f64&gt;::zeros(lworkl);
    let mut info = 0;

    let sigma = 0.0;
    loop
    {
        unsafe {
            dsaupd_c(
                &amp;mut ido,
                &quot;I&quot;.as_ptr() as *const i8,
                dimension as i32,
                eigenvalue_kind.as_ptr() as *const i8,
                number_of_eigenvalues as i32,
                f64::EPSILON,
                residual.as_ptr_mut() as *mut f64,
                number_of_basis_vectors as i32,
                arnoldi_vectors.as_ptr_mut() as *mut f64,
                dimension as i32,
                iparam.as_mut_ptr(), // no idea what this does
                ipntr.as_mut_ptr(),  // no idea what this does
                workd.as_ptr_mut() as *mut f64,
                workl.as_ptr_mut() as *mut f64,
                lworkl as i32,
                &amp;mut info,
            );

            match info
            {
                0 | 1 | 2 =&gt;
                {}
                -1 =&gt; panic!(&quot;N must be positive.&quot;),
                -2 =&gt; panic!(&quot;NEV must be positive.&quot;),
                -3 =&gt; panic!(&quot;NCV-NEV &gt;= 2 and less than or equal to N.&quot;),
                -4 =&gt; panic!(&quot;Maximum iterations must be greater than 0.&quot;),
                -5 =&gt; panic!(&quot;Maximum iterations must be greater than 0.&quot;),
                i =&gt; panic!(&quot;dsaupd_c returned error code {}&quot;, i),
            }

            if (ido == -1) || (ido == 1)
            {
                let res = mat * workd.subrows(ipntr[0] as usize - 1, dimension);
                let y = faer_shift_invert_solve(mat, sigma, &amp;res);
                workd
                    .subrows_mut(ipntr[1] as usize - 1, dimension)
                    .copy_from(&amp;y);
            }
            else if ido == 99
            {
                break;
            }
            else
            {
                panic!(&quot;ido code: {}&quot;, ido);
            }
        }
    }

    let select = vec![false as i32; number_of_basis_vectors];
    let mut d: Col&lt;f64&gt; = Col::zeros(number_of_eigenvalues + 1);
    let mut z: Mat&lt;f64&gt; = Mat::zeros(dimension, number_of_basis_vectors);
    unsafe {
        dseupd_c(
            vectors as i32,
            &quot;A&quot;.as_ptr() as *const i8,
            select.as_ptr(),
            d.as_ptr_mut() as *mut f64,
            z.as_ptr_mut() as *mut f64,
            dimension as i32,
            0.0 as f64,
            &quot;I&quot;.as_ptr() as *const i8,
            dimension as i32,
            eigenvalue_kind.as_ptr() as *const i8,
            number_of_eigenvalues as i32,
            f64::EPSILON,
            residual.as_ptr_mut() as *mut f64,
            number_of_basis_vectors as i32,
            arnoldi_vectors.as_ptr_mut() as *mut f64,
            dimension as i32,
            iparam.as_mut_ptr(),
            ipntr.as_mut_ptr(),
            workd.as_ptr_mut() as *mut f64,
            workl.as_ptr_mut() as *mut f64,
            lworkl as i32,
            &amp;mut info,
        );
    }

    match info {
        0 =&gt; { /* Normal exit, do nothing */ }
        1 =&gt; panic!(&quot;INFO = 1: The Schur form computed by LAPACK routine dsytrd failed to reorder or to compute eigenvalues and eigenvectors to high relative accuracy.&quot;),
        3 =&gt; panic!(&quot;INFO = 3: No shifts could be applied during a cycle of the implicitly restarted Arnoldi iteration. One possible reason is that the Schur vectors failed to converge.&quot;),
        -1 =&gt; panic!(&quot;INFO = -1: N must be greater than 0.&quot;),
        -2 =&gt; panic!(&quot;INFO = -2: NEV must be greater than 0 and less than N.&quot;),
        -3 =&gt; panic!(&quot;INFO = -3: NCV must be greater than NEV and less than or equal to N.&quot;),
        -5 =&gt; panic!(&quot;INFO = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'.&quot;),
        -6 =&gt; panic!(&quot;INFO = -6: BMAT must be one of 'I' or 'G'.&quot;),
        -7 =&gt; panic!(&quot;INFO = -7: The part of the spectrum specified by WHICH is not supported.&quot;),
        -8 =&gt; panic!(&quot;INFO = -8: The spectrum range parameters are inconsistent.&quot;),
        -9 =&gt; panic!(&quot;INFO = -9: WORKD[1:N] is not correctly allocated.&quot;),
        -10 =&gt; panic!(&quot;INFO = -10: The starting vector is zero.&quot;),
        -11 =&gt; panic!(&quot;INFO = -11: The standard Lanczos iteration (for generalized symmetric eigenvalue problem) failed to converge.&quot;),
        -12 =&gt; panic!(&quot;INFO = -12: The shifts could not be applied during a cycle of the implicitly restarted Arnoldi iteration. One possible reason is that the Schur vectors failed to converge.&quot;),
        _ =&gt; panic!(&quot;Unknown INFO value: {}&quot;, info),
    }

    let mut eigs = d.subrows(0, number_of_eigenvalues).to_owned();
    for eig in eigs.iter_mut()
    {
        *eig = sigma + 1.0 / *eig;
    }

    (eigs, z)
}

fn faer_shift_invert_solve(
    a: &amp;SparseColMat&lt;usize, f64&gt;,
    sigma: f64,
    b: &amp;Col&lt;f64&gt;,
) -&gt; Col&lt;f64&gt;
{
    let triplets = (0..a.nrows()).map(|i| (i, i, 1.)).collect::&lt;Vec&lt;_&gt;&gt;();
    let identity =
        SparseColMat::try_new_from_triplets(a.nrows(), a.nrows(), &amp;triplets).unwrap();

    use faer::prelude::SpSolver;
    let qr = (a - sigma * identity).sp_qr().unwrap();
    qr.solve(b)
}
</code></pre>

CC BY-SA 4.0
2024-07-22T00:08:28.517Z
Potentially related: https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/87
https://docs.scipy.org/doc/scipy/tutorial/arpack.html
If you are using the cotan Laplacian you know the zero eigenvalue has corresponding constant eigenvectors, so you could just avoid computing it with arpack. Another option is to add $\frac{s}{n}11^T$ to your matrix, then it will be non-singular with eigenvalue $s$ now for the constant eigenvector, you could make $s$ large enough as to move it away from the small ones so arpack can compute the relevant ones. I have no ideas about geogram.
lightxbulb
CC BY-SA 4.0
2024-07-22T00:19:08.430Z
Use the sift-invert option instead of SM as described in the scipy tutorial. You can potentially couple this with the modification of the matrix so you also shift out the zero eigenvalue. And it would be nice if you could post some pretty pictures of a visualization of the eigenfunctions if that works out.
lightxbulb
CC BY-SA 4.0
2024-07-22T05:17:53.173Z
@lightxbulb So I am reading both the scipy and arpack's documentation, and i am finding a problem. I am working on rust, so I don;t have access to eigsh. And eigsh seems to automagically define the M matrix needed for the invert and shift trick.

Do you know what I should do to construct that matrix to apply the invert shift operation?
Makogan
CC BY-SA 4.0
2024-07-22T08:17:49.550Z
Both eigs and eigsh have a sigma parameter which is what you need for the shift. You don't need to specify $M$ since you are not solving a generalized eigenvalue problem (then $M=I$ is used). Arpack will solve $(A-\sigma I)^{-1}x = \nu x$ for the largest magnitude eigs, i.e., $\frac{1+\sigma \nu}{\nu}x = Ax$, $\lambda = \frac{1}{\nu} +\sigma$, $\nu =\frac{1}{\lambda - \sigma}$. These are the smallest eigs of $A$ around $\sigma$. That said using eigs is probably highly suboptimal for symmetric matrices as it likely uses Arnoldi, unlike eigsh which should use Lanczos.
lightxbulb
CC BY-SA 4.0
2024-07-22T17:04:04.503Z
@lightxbulb I think we might be misscomunicating. eigs and eigsh are not part of arpack. They are higher level utilities provided byt the python package your provided. arpack only has *aupd and *eupd subroutines.
Makogan
CC BY-SA 4.0
2024-07-22T17:27:39.163Z
Star * in aupd and eupd is {s,d,c,z}{s,n}, s = single precision float, d = double precision float, c = single precision complex, z = double precision comex. The second: s = symmetric/Hermitian, n = non-symmetric. So eigsh is referring so the second flag being s, while eigs is referring to second flag being n. Python calls these under the hood. See table A.2. here for details: http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf
If you want to use shift-invert, you should provide a routine that can solve $(A-\sigma I)x = y$.
lightxbulb
CC BY-SA 4.0
2024-07-22T17:42:37.530Z
@lightxbulb I am using shift-invert and I am still getting the ame weird eigenvectors. Worse, to try to test things out I made a sparse symmetric tri-fiagonal toeplitz matrix, so that I knew what the eigenvalues are a priori and I am getting nothing remotely close to the expected eigenvalues.
Makogan
CC BY-SA 4.0
2024-07-22T18:15:21.010Z
You can try this on e.g. a 10x10 Laplacian matrix, and compare to e.g. matlabs eigs or scipy's eigs, if it's still wrong I suppose you have a bug in your implementation. I would assume the bug is not in arpack, although maybe there could be some issues with singular matrices.
lightxbulb
CC BY-SA 4.0
2024-07-23T20:30:27.883Z
I did as suggested. My eigenvalues are entirely different from what scipy returns so there is something wrong. But I am unsure where. It has to be in the communication between arpack and the surrounding rust code as I trust the library bring used for the linear algebra operations. But I don;t know if it is a wrong pointer, data miss alignment or anything else of that ilk.
Makogan
<computational-geometry><geometry><matrix><linear-algebra>
Manifold harmonics vectors not working?
Makogan
# Manifold harmonics vectors not working?

(Asked by: username_0 on 2024-07-21)

I am trying to implement [manifold harmonics](https://graphics.stanford.edu/courses/cs233-21-spring/ReferencedPapers/understand_geometry_01631196.pdf).

To see if I am doing things right I am visualizing the [Fiedler eigenvector](https://en.wikipedia.org/wiki/Algebraic_connectivity), which is the eigenvector associated with the second smallest eigenvalue.

On my own code base (which is using ARPACK to compute the eigendecomposition), I am getting strange results.

Normally, if I am requesting the k smallest eigenvalues, according to the wikipedia page for a graph laplacian, the first eigenvalue should be 0. However, in my case I am getting these eigenvalues:

```rs
[-4211.039032568338],
[2364.2174825349684],
[3668.6929511652374],
[11190.921998283911],
[17989.923028566314],

```

The visualization of the eigenvectors is also weird for example I am getting this image for the second eigenvector (sorted from smallest eigenvalue to largest).

[![enter image description here](https://i.sstatic.net/TMbDbleJ.png)](https://i.sstatic.net/TMbDbleJ.png)

I am not sure if it is the input matrix or my use of arpack that is breaking, I would be more inclined of thinking it is arpack because i have tested the matrix in the past.

To try to debug I build geogram to visualize their results for the same mesh. But it doesn't seem to do anything useful.

I loaded the mesh:

[![enter image description here](https://i.sstatic.net/LcXILBdr.png)](https://i.sstatic.net/LcXILBdr.png)

Clicked on manifold harmonics:

[![enter image description here](https://i.sstatic.net/4r5p8TLj.png)](https://i.sstatic.net/4r5p8TLj.png)

Applied them:

[![enter image description here](https://i.sstatic.net/w6mI16Y8.png)](https://i.sstatic.net/w6mI16Y8.png)

For every attribute labeled MH I am geting that the range is 0

[![enter image description here](https://i.sstatic.net/AcVKrx8J.png)](https://i.sstatic.net/AcVKrx8J.png)

As if the eigenvectors are always 0. I have tried multiple meshes and they all return this result. I was hoping someone could help me.

As suggested in the comments, I am attempting to use shift-invert:

```rs
pub fn eigen_decompose_sparse_matrix(
    mat: &SparseColMat<usize, f64>,
) -> (Col<f64>, Mat<f64>)
{
    debug_assert!(mat.nrows() == mat.ncols());
    faer_arpack_symmetric_f64(
        &mat,
        mat.nrows(),
        Eigenvaluekind::LargestMagnitude.as_str(),
        5,
        10,
        3000,
        true,
    )
}

pub enum Eigenvaluekind
{
    LargestMagnitude,
    SmallestMagnitude,
    LargestRealPart,
    SmallestRealPart,
    LargestImaginaryPart,
    SmallestImaginaryPart,
}

impl Eigenvaluekind
{
    fn as_str(&self) -> &'static str
    {
        match self
        {
            Eigenvaluekind::LargestMagnitude => "LM",
            Eigenvaluekind::SmallestMagnitude => "SM",
            Eigenvaluekind::LargestRealPart => "LR",
            Eigenvaluekind::SmallestRealPart => "SR",
            Eigenvaluekind::LargestImaginaryPart => "LI",
            Eigenvaluekind::SmallestImaginaryPart => "SI",
        }
    }
}

// http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf
fn faer_arpack_symmetric_f64(
    mut mat: &SparseColMat<usize, f64>,
    dimension: usize,
    eigenvalue_kind: &str,
    number_of_eigenvalues: usize,
    number_of_basis_vectors: usize,
    maxiter: usize,
    vectors: bool,
) -> (Col<f64>, Mat<f64>)
{
    let mut ido = 0;
    let mut residual: Col<f64> = Col::zeros(dimension);
    let mut arnoldi_vectors = Mat::<f64>::zeros(dimension, number_of_basis_vectors);
    let mut iparam = [0; 11];
    iparam[0] = 1;
    iparam[2] = maxiter as i32;
    iparam[6] = 1;
    let mut ipntr = [0; 14];
    let mut workd = Col::zeros(3 * dimension);
    let lworkl = number_of_basis_vectors * (number_of_basis_vectors + 8);
    let mut workl = Col::<f64>::zeros(lworkl);
    let mut info = 0;

    let sigma = 0.0;
    loop
    {
        unsafe {
            dsaupd_c(
                &mut ido,
                "I".as_ptr() as *const i8,
                dimension as i32,
                eigenvalue_kind.as_ptr() as *const i8,
                number_of_eigenvalues as i32,
                f64::EPSILON,
                residual.as_ptr_mut() as *mut f64,
                number_of_basis_vectors as i32,
                arnoldi_vectors.as_ptr_mut() as *mut f64,
                dimension as i32,
                iparam.as_mut_ptr(), // no idea what this does
                ipntr.as_mut_ptr(),  // no idea what this does
                workd.as_ptr_mut() as *mut f64,
                workl.as_ptr_mut() as *mut f64,
                lworkl as i32,
                &mut info,
            );

            match info
            {
                0 | 1 | 2 =>
                {}
                -1 => panic!("N must be positive."),
                -2 => panic!("NEV must be positive."),
                -3 => panic!("NCV-NEV >= 2 and less than or equal to N."),
                -4 => panic!("Maximum iterations must be greater than 0."),
                -5 => panic!("Maximum iterations must be greater than 0."),
                i => panic!("dsaupd_c returned error code {}", i),
            }

            if (ido == -1) || (ido == 1)
            {
                let res = mat * workd.subrows(ipntr[0] as usize - 1, dimension);
                let y = faer_shift_invert_solve(mat, sigma, &res);
                workd
                    .subrows_mut(ipntr[1] as usize - 1, dimension)
                    .copy_from(&y);
            }
            else if ido == 99
            {
                break;
            }
            else
            {
                panic!("ido code: {}", ido);
            }
        }
    }

    let select = vec![false as i32; number_of_basis_vectors];
    let mut d: Col<f64> = Col::zeros(number_of_eigenvalues + 1);
    let mut z: Mat<f64> = Mat::zeros(dimension, number_of_basis_vectors);
    unsafe {
        dseupd_c(
            vectors as i32,
            "A".as_ptr() as *const i8,
            select.as_ptr(),
            d.as_ptr_mut() as *mut f64,
            z.as_ptr_mut() as *mut f64,
            dimension as i32,
            0.0 as f64,
            "I".as_ptr() as *const i8,
            dimension as i32,
            eigenvalue_kind.as_ptr() as *const i8,
            number_of_eigenvalues as i32,
            f64::EPSILON,
            residual.as_ptr_mut() as *mut f64,
            number_of_basis_vectors as i32,
            arnoldi_vectors.as_ptr_mut() as *mut f64,
            dimension as i32,
            iparam.as_mut_ptr(),
            ipntr.as_mut_ptr(),
            workd.as_ptr_mut() as *mut f64,
            workl.as_ptr_mut() as *mut f64,
            lworkl as i32,
            &mut info,
        );
    }

    match info {
        0 => { /* Normal exit, do nothing */ }
        1 => panic!("INFO = 1: The Schur form computed by LAPACK routine dsytrd failed to reorder or to compute eigenvalues and eigenvectors to high relative accuracy."),
        3 => panic!("INFO = 3: No shifts could be applied during a cycle of the implicitly restarted Arnoldi iteration. One possible reason is that the Schur vectors failed to converge."),
        -1 => panic!("INFO = -1: N must be greater than 0."),
        -2 => panic!("INFO = -2: NEV must be greater than 0 and less than N."),
        -3 => panic!("INFO = -3: NCV must be greater than NEV and less than or equal to N."),
        -5 => panic!("INFO = -5: WHICH must be one of 'LM', 'SM', 'LR', 'SR', 'LI', 'SI'."),
        -6 => panic!("INFO = -6: BMAT must be one of 'I' or 'G'."),
        -7 => panic!("INFO = -7: The part of the spectrum specified by WHICH is not supported."),
        -8 => panic!("INFO = -8: The spectrum range parameters are inconsistent."),
        -9 => panic!("INFO = -9: WORKD[1:N] is not correctly allocated."),
        -10 => panic!("INFO = -10: The starting vector is zero."),
        -11 => panic!("INFO = -11: The standard Lanczos iteration (for generalized symmetric eigenvalue problem) failed to converge."),
        -12 => panic!("INFO = -12: The shifts could not be applied during a cycle of the implicitly restarted Arnoldi iteration. One possible reason is that the Schur vectors failed to converge."),
        _ => panic!("Unknown INFO value: {}", info),
    }

    let mut eigs = d.subrows(0, number_of_eigenvalues).to_owned();
    for eig in eigs.iter_mut()
    {
        *eig = sigma + 1.0 / *eig;
    }

    (eigs, z)
}

fn faer_shift_invert_solve(
    a: &SparseColMat<usize, f64>,
    sigma: f64,
    b: &Col<f64>,
) -> Col<f64>
{
    let triplets = (0..a.nrows()).map(|i| (i, i, 1.)).collect::<Vec<_>>();
    let identity =
        SparseColMat::try_new_from_triplets(a.nrows(), a.nrows(), &triplets).unwrap();

    use faer::prelude::SpSolver;
    let qr = (a - sigma * identity).sp_qr().unwrap();
    qr.solve(b)
}

```

**Tags:** computational-geometry, geometry, matrix, linear-algebra

> **Comments:**
> **username_1**: Potentially related: https://github.com/JuliaLinearAlgebra/Arpack.jl/issues/87 https://docs.scipy.org/doc/scipy/tutorial/arpack.html If you are using the cotan Laplacian you know the zero eigenvalue has corresponding constant eigenvectors, so you could just avoid computing it with arpack. Another option is to add $\frac{s}{n}11^T$ to your matrix, then it will be non-singular with eigenvalue $s$ now for the constant eigenvector, you could make $s$ large enough as to move it away from the small ones so arpack can compute the relevant ones. I have no ideas about geogram.  
> **username_1**: Use the sift-invert option instead of SM as described in the scipy tutorial. You can potentially couple this with the modification of the matrix so you also shift out the zero eigenvalue. And it would be nice if you could post some pretty pictures of a visualization of the eigenfunctions if that works out.  
> **username_0**: @username_1 So I am reading both the scipy and arpack's documentation, and i am finding a problem. I am working on rust, so I don;t have access to eigsh. And eigsh seems to automagically define the M matrix needed for the invert and shift trick.  Do you know what I should do to construct that matrix to apply the invert shift operation?  
> **username_1**: Both eigs and eigsh have a sigma parameter which is what you need for the shift. You don't need to specify $M$ since you are not solving a generalized eigenvalue problem (then $M=I$ is used). Arpack will solve $(A-\sigma I)^{-1}x = \nu x$ for the largest magnitude eigs, i.e., $\frac{1+\sigma \nu}{\nu}x = Ax$, $\lambda = \frac{1}{\nu} +\sigma$, $\nu =\frac{1}{\lambda - \sigma}$. These are the smallest eigs of $A$ around $\sigma$. That said using eigs is probably highly suboptimal for symmetric matrices as it likely uses Arnoldi, unlike eigsh which should use Lanczos.  
> **username_0**: @username_1 I think we might be misscomunicating. eigs and eigsh are not part of arpack. They are higher level utilities provided byt the python package your provided. arpack only has *aupd and *eupd subroutines.  
> **username_1**: Star * in aupd and eupd is {s,d,c,z}{s,n}, s = single precision float, d = double precision float, c = single precision complex, z = double precision comex. The second: s = symmetric/Hermitian, n = non-symmetric. So eigsh is referring so the second flag being s, while eigs is referring to second flag being n. Python calls these under the hood. See table A.2. here for details: http://li.mit.edu/Archive/Activities/Archive/CourseWork/Ju_Li/MITCourses/18.335/Doc/ARPACK/Lehoucq97.pdf If you want to use shift-invert, you should provide a routine that can solve $(A-\sigma I)x = y$.  
> **username_0**: @username_1 I am using shift-invert and I am still getting the ame weird eigenvectors. Worse, to try to test things out I made a sparse symmetric tri-fiagonal toeplitz matrix, so that I knew what the eigenvalues are a priori and I am getting nothing remotely close to the expected eigenvalues.  
> **username_1**: You can try this on e.g. a 10x10 Laplacian matrix, and compare to e.g. matlabs eigs or scipy's eigs, if it's still wrong I suppose you have a bug in your implementation. I would assume the bug is not in arpack, although maybe there could be some issues with singular matrices.  
> **username_0**: I did as suggested. My eigenvalues are entirely different from what scipy returns so there is something wrong. But I am unsure where. It has to be in the communication between arpack and the surrounding rust code as I trust the library bring used for the linear algebra operations. But I don;t know if it is a wrong pointer, data miss alignment or anything else of that ilk.  

---

## 1 answers

### Answer by username_0 on 2024-07-26. Score: 1 | Accepted answer

So there were multiple issues with my use of ARPACK. Some of them:

Linear algebra structs not having continous memory alignment, improper use of the shift invert mode, not having enough lanczos vectors...

This is a result after having fixed them:

[![enter image description here](https://i.sstatic.net/pzK79tAf.png)](https://i.sstatic.net/pzK79tAf.png)

In more detail:

Faer's matrices are not guaranteed to be stored in contiguous memory, so using them for high level storage of the results from ARPACK required some care to avoid unaligning the memory.

For the shift invert mode, I was not setting the correct code to the iparam routine.

And I had neglected to read the part of the documentation that suggests having about 3x more lanczos basis vectors than requested eigenvalues.

> **Comments:**
> **username_1**: It may be useful to others working with rust if you provide more details on how you solved the three problems you mentioned.  
> **username_0**: @username_1 added  

---

