# Get an eigvec by eigval, using eigvecs, only allow SymTridiagonal?

**URL:** https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408
**Category:** New to Julia
**Tags:** question, linearalgebra
**Created:** [March 27, 2025, 11:14am UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408 "2025-03-27T11:14:15Z")
**Posts on this page:** 13
**Page:** 1

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [March 27, 2025, 11:14am UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/1 "2025-03-27T11:14:15Z")

</div>

Hi, there

I wonder why the `eigvecs` method does not work for a `A::Symmetric{Int64, Matrix{Int64}}`.

```julia
using LinearAlgebra
A = SymTridiagonal([1.; 2.; 1.], [2.; 3.])
val3 = eigvals(A)[3]
vec3 = reshape(eigvecs(A, [val3]), (:,)) # ✅ this works
A * vec3 .- val3 * vec3 # double check the correctness
# re-do the above prosess
A = Symmetric(Matrix(A))

julia> eigvecs(A, [val3])
ERROR: MethodError: no method matching eigvecs(::Symmetric{Int64, Matrix{Int64}}, ::Vector{Float64})
The function `eigvecs` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  eigvecs(::AbstractMatrix, ::AbstractMatrix; kws...)
   @ LinearAlgebra K:\julia-1.11.4\share\julia\stdlib\v1.11\LinearAlgebra\src\eigen.jl:651
  eigvecs(::SymTridiagonal{<:Union{Float32, Float64, ComplexF64, ComplexF32}, <:StridedVector{T} where T}, ::Vector{<:Real})
   @ LinearAlgebra K:\julia-1.11.4\share\julia\stdlib\v1.11\LinearAlgebra\src\tridiag.jl:326
  eigvecs(::Union{Hermitian{T, S}, Symmetric{T, S}} where {T, S})
   @ LinearAlgebra K:\julia-1.11.4\share\julia\stdlib\v1.11\LinearAlgebra\src\symmetriceigen.jl:251
  ...

```

Should I use the following method instead? Is this the correct usage?

```julia
function my_eigvecs(A::Symmetric, val) return nullspace(A - val * one(A))[:, 1] end
vec3 = my_eigvecs(A, val3) # retrieve an eigenvector

```

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 27, 2025, 1:12pm UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/2 "2025-03-27T13:12:07Z")

</div>

The `eigvecs(A, λ)` method is currently only implemented for `SymTridiagonal` matrices. For other matrix types there is currently only `eigvecs(A)`.

In principle, we could easily extend this to real-symmetric/Hermitian matrices, since they can be transformed into real `SymTridiagonal` matrices by the Hessenberg factorization. This should work:

```julia
using LinearAlgebra
import LinearAlgebra: eigvecs, RealHermSymComplexHerm

function eigvecs(A::RealHermSymComplexHerm, λ::AbstractVector{<:Real})
    F = hessenberg(A) # transform to SymTridiagonal form
    X = eigvecs(F.H, λ)
    return F.Q * X # transform eigvecs of F.H back to eigvecs of A
end

```

For only computing a small subset of the eigenvectors, it seems to be a couple times faster than `eigvecs`, especially if you don’t count the cost of the eigenvalues (e.g. you already have them for some other reason).

Might be worth putting together a PR to [LinearAlgebra.jl](https://github.com/JuliaLang/LinearAlgebra.jl) if you are interested in this functionality? [eigvecs(A::Hermitian, eigvals) method? · Issue #1248 · JuliaLang/LinearAlgebra.jl · GitHub](https://github.com/JuliaLang/LinearAlgebra.jl/issues/1248)

> [@WalterMadelim](#):
>
> Should I use the following method instead? Is this the correct usage?

No. `nullspace` employs an SVD, which is as costly as computing all the eigenvectors and eigenvalues with `eigen`. (In fact, the SVD calls `eigen` for Hermitian matrices.)

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 27, 2025, 1:39pm UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/3 "2025-03-27T13:39:04Z")

</div>

> [@stevengj](#):
>
> In principle, we could easily extend this to real-symmetric/Hermitian matrices, since they can be transformed into real `SymTridiagonal` matrices by the Hessenberg factorization.

As [I commented in the LinearAlgebra.jl issue](https://github.com/JuliaLang/LinearAlgebra.jl/issues/1248#issuecomment-2758065757), however, if you are calling `eigvals(A)` explicitly, you are better off doing the Hessenberg factorization yourself and re-using it for _both_ `eigvecs` and `eigvals`, since the Hessenberg factorization is the most expensive part of Hermitian eigensolves.

(I really should write a blog post about Hessenberg factorizations in Julia at some point; most people aren’t familiar with this factorization and why it is important.)

---

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [March 28, 2025, 1:07am UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/4 "2025-03-28T01:07:32Z")

</div>

The motivation for _directly_ retrieving a eigenvector associated with the minimum eigenvalue is here:

Assume `A` is a symmetric real matrix whose `size(A) == (n, n)`.  
According to the definition, if `A` is _not_ PSD, then we should could find a `v`, such that

```julia
transpose(v) * A * v == eigmin(A) < 0

```

Once we know `v`, we could add the linear cut `@constraint(model, transpose(v) * X * v >= 0)`, as one member of the full PSD cut `@constraint(model, X in PSDCone())`.

Therefore, I hope that `LinearAlgebra` could extend the `eigmin` function such that it **also provide a corresponding eigenvector**.

---

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [March 28, 2025, 1:29am UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/5 "2025-03-28T01:29:32Z")

</div>

> [@stevengj](#):
>
> if you are interested in this functionality?

Yes, it is useful, as my motivation above.

> [@stevengj](#):
>
> most people aren’t familiar with this factorization

Ingenuous students probably didn’t learn them well at school. If they are, they should spend more effort learning by themselves. If they learn pure knowledge without applications, they may forget all of them after e.g. 1 year. 🙂

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 28, 2025, 3:40am UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/6 "2025-03-28T03:40:29Z")

</div>

> [@WalterMadelim](#):
>
> Therefore, I hope that `LinearAlgebra` could extend the `eigmin` function such that it **also provide a corresponding eigenvector**.

You can just call `eigen(A, 1:1)`, no?

---

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [March 28, 2025, 4:06am UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/7 "2025-03-28T04:06:58Z")

</div>

Ah, yes. But is it efficient? x-ref [this](https://github.com/JuliaLang/LinearAlgebra.jl/issues/1248#issuecomment-2760008073).

---

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [March 28, 2025, 4:16am UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/8 "2025-03-28T04:16:12Z")

</div>

> [@stevengj](#):
>
> eigen(A, 1:1)

This API works, but is somewhat not intelligible for the user.  
Why can’t it be `eigen(A, begin:begin)`, or if I want to retrieve the maximum,  
`eigen(A, end:end)`. These cannot work out.

Now I’m updating like this

```julia
import LinearAlgebra
function my_eigmin(A)
    valmin, vecmin = LinearAlgebra.eigen(A, 1:1)
    return first(valmin), vecmin[:, 1]
end
# test code begins
A = LinearAlgebra.SymTridiagonal([1.; 2.; 1.], [2.; 3.])
A = LinearAlgebra.Symmetric(Matrix(A))
valmin, vecmin = my_eigmin(A)

```

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [March 28, 2025, 12:49pm UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/9 "2025-03-28T12:49:39Z")

</div>

> [@WalterMadelim](#):
>
> Ah, yes. But is it efficient? x-ref [this](https://github.com/JuliaLang/LinearAlgebra.jl/issues/1248#issuecomment-2760008073).

Yes. `eigen(A, 1:1)` only works for Hermitian matrices, where there is a specialized algorithm. The method you are referring to in that comment is `eigmin` for general matrices, which is not as efficient.

---

<div class="post-metadata">

### Author: ![jishnub](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jishnub/32/33620_2.png) [@jishnub](https://discourse.julialang.org/u/jishnub)
#### Post date: [March 28, 2025, 1:08pm UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/10 "2025-03-28T13:08:40Z")

</div>

> [@WalterMadelim](#):
>
> Why can’t it be `eigen(A, begin:begin)`, or if I want to retrieve the maximum,  
> `eigen(A, end:end)`.

`begin` and `end` are allowed to point to the first and last indices only within an indexing operation. There simply isn’t a way to pass these generically as arguments. It might be possible to avail a somewhat similar syntax by using EndpointRanges.jl.

---

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [May 1, 2025, 2:53pm UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/11 "2025-05-01T14:53:48Z")

</div>

> [@stevengj](#):
>
> about Hessenberg factorizations in Julia

I find a noteworthy point, written in the comment below

```julia
using LinearAlgebra

N = rand(2:17) 
A = rand(-9:0.17:9, N, N)
F = hessenberg(A)
Q = Matrix(F.Q) # the first column of Q is (artificially) fixed?
@assert Q[1, :] == Q[:, 1] == Matrix(I, N, N)[:, 1]

```

Is this generally true? (The algorithm is from LAPACK?)  
Does this setting makes the (Hessenberg reduction) unique?  
**Edit** : well, if a sequence of Householder matrices is employed, then each component Householder matrix is bound to has this characteristic, and so has their product (i.e. the `Q`)

Nonetheless, for a `Symmetric` matrix `A`, the situation is the opposite

```julia
using LinearAlgebra

N = rand(2:17) 
A = Symmetric(rand(-9:0.17:9, N, N))
F = hessenberg(A)
Q = Matrix(F.Q)
@assert Q[end, :] == Q[:, end] == Matrix(I, N, N)[:, end]

```

This might due to the action of zeroing elements from bottom to top.

---

<div class="post-metadata">

### Author: ![WalterMadelim](https://avatars.discourse-cdn.com/v4/letter/w/3e96dc/32.png) [@WalterMadelim](https://discourse.julialang.org/u/WalterMadelim)
#### Post date: [May 3, 2025, 1:05pm UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/12 "2025-05-03T13:05:59Z")

</div>

About retrieving the smallest eigenvalue along with an eigenvector associated with it, I find this code in LinearAlgebra.jl

```julia
eigvals!(A::SymTridiagonal{<:BlasReal,<:StridedVector}, irange::UnitRange) =
    LAPACK.stegr!('N', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)[1]

```

I don’t know what is the underlying algorithm. But I have a rough guess—using power iteration, then inverse iteration.

Assume `A` is a real symmetric matrix, then it has `n` real eigenvalues.  
We can firstly use power iteration to derive its dominant eigenvalue.

1. If the outcome happens to be negative, then we quit.
2. Otherwise, at this time we know the modulus of the dominant eigenvalue. The negative counterpart of this value can serve as an initial guess of the smallest eigenvalue. Then we can perform inverse iteration. After that, we quit.

A simple code associated with the idea above is as follows

```julia
using LinearAlgebra
function get_initial_q()
    q = rand(-1:.0017:1, 3)
    q /= norm(q)
    q
end
macro iteration_code()
    esc(quote
        q = get_initial_q()
        for i in 1:50
            t = A * q
            q = t / norm(t)
            @info "[$i]: monitor" q
            sleep(.1)
        end
    end)
end
# we are interested in the `-3` below, but we don't know it ab initio
A = Diagonal([1., -3, 4])
@iteration_code() # [Power Iteration] this reveals the modulus of the dominant eigenvalue, which is `4`
# we set `-4` as the initial guess of the smallest eigenvalue  
A = inv(A - (-4) * I) # only shows the idea, not applicable in practice
@iteration_code() # [Inverse Iteration]

```

If there is anyone who have more knowledge on this, please correct me. 🙂

---

<div class="post-metadata">

### Author: ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 29, 2025, 10:40pm UTC](https://discourse.julialang.org/t/get-an-eigvec-by-eigval-using-eigvecs-only-allow-symtridiagonal/127408/13 "2025-05-29T22:40:55Z")

</div>

> [@WalterMadelim](#):
>
> I don’t know what is the underlying algorithm. But I have a rough guess—using power iteration, then inverse iteration.

No, it looks like it calls [`dstemr`](https://netlib.org/lapack/explore-html/d4/dec/group__stemr_ga71cbf49a0387762afa39582e6abbe466.html#ga71cbf49a0387762afa39582e6abbe466) under the hood, which says:

> Depending on the number of desired eigenvalues, these are computed either by bisection or the dqds algorithm.

> [@WalterMadelim](#):
>
> We can firstly use power iteration to derive its dominant eigenvalue.

Power iterations converge slowly if you have nearby eigenvalues; you should really only use them if you have some knowledge of the spectrum.
