Taking Eigen seriously #48549

antoine-levitt · 2023-02-06T10:29:15Z

This extends the Eigen API to support a few more functions (matvecs, matrix functions), and optimizes in the case where the eigenvectors are unitary. This is very convenient for eg computing exp(tA)x at different times.

dkarrasch · 2023-02-06T14:22:29Z

x-ref #41830

antoine-levitt · 2023-02-06T14:34:44Z

Oh wow, I see I even commented on that PR. It's not even that long ago, and I completely forgot about it... I'd say the present PR is simpler, since it does not introduce extra API (new functions or arguments) but only add dispatches for things that "should" work. It is also safer (but slightly more restrictive) because it doesn't allow eigen to set normality, but only sets it when called on hermitian matrices (the most common case)

dkarrasch · 2023-02-06T15:23:13Z

but only sets it when called on hermitian matrices (the most common case)

I believe we wrap Hermitian matrices with Hermitian when we detect it? So that should be a good enough option.

I just wanted to raise awareness for that other PR, so that "overlapping" aspects are well-resolved. Would be awesome to also push the other PR across the finish line. @RalphAS

antoine-levitt · 2023-02-07T08:38:42Z

I just wanted to raise awareness for that other PR, so that "overlapping" aspects are well-resolved. Would be awesome to also push the other PR across the finish line. @RalphAS

My PR extends the Eigen type by adding dispatches for matrix functions, matrix multiplications, etc and only adds a "normal" field to the Eigen type. @RalphAS 's adds the "normal" field, adds a new matrix_function and adds more optional computations (left vectors, condition numbers) to eigen (stored in the Eigen type). So really, apart from the addition of the normal field there's not that much overlap, and upon merging one the other can be adapted easily. I would argue in favor of my approach to matrix functions as being more idiomatic, though. It doesn't support arbitrary functions, but for that I think we really need a MatrixFunctions.jl package that supports non-symmetric matrices anyway.

RalphAS · 2023-02-08T04:45:56Z

I like this PR so far. This is a good opportunity to document some of the "hidden" capabilities of the Eigen struct which you have extended.

I think the bulk of #41830 really needs new data structures that haven't been worked out yet, so it should be orthogonal.

antoine-levitt · 2023-02-09T08:25:39Z

This is a good opportunity to document some of the "hidden" capabilities of the Eigen struct which you have extended.

Good point. I added news and docs.

dkarrasch

LGTM. Shall we rename are_vectors_unitary to unitary_[e]vecs or something?

antoine-levitt · 2023-02-09T14:49:16Z

I don't know, it's a pretty internal field so I feel being super explicit is good and long identifiers are not an issue, unitary_evecs might make somebody think it's an array of vectors that are more unitary than vectors or something. But I'm fine to rename to whatever if you want.

garrison · 2023-02-12T07:55:16Z

NEWS.md

@@ -47,6 +47,9 @@ Standard library changes
  supports precompilation and setting compatibility.
 #### LinearAlgebra

+- The `Eigen` type now supports more methods, such as matrix functions
+  (which return `Eigen` objects), adjoints and matrix-vector product.
+  `inv(::Eigen)` now produces an `Eigen` type instead of a matrix.


Is this enough for this PR to gain the "minor change" label?

FWIW, I did an audit of my personal code, and while I never assume that inv(::Eigen) returns a matrix, I have relied on the fact that inv(::LU) returns a matrix.

In any case, it might be nice to say explicitly how one can continue to get the inverse matrix from ::Eigen, if that's what a user wants.

Good point. x-ref #32137. And certainly this requires a pkgeval run if we go with that change.

As far as I understood the rules were more relaxed about breaking changes for stdlibs? I can either special-case inv to return a matrix (and we can change it for 2.0) or make the change now with a more explicit NEWS. Can somebody run pkgeval, just to get an idea of how many people rely on this? (I'd be surprised, but...)

@antoine-levitt You can run pkgeval yourself! It's as easy as:

@nanosoldier runtests()

See https://github.com/JuliaCI/Nanosoldier.jl#readme for more info

dkarrasch · 2023-02-12T10:15:52Z

Another issue is that eigenvalues are no longer ordered in ascending order, as we typically assume for values, certainly for symmetric/Hermitian matrices...

antoine-levitt · 2023-02-13T08:02:45Z

Another issue is that eigenvalues are no longer ordered in ascending order, as we typically assume for values, certainly for symmetric/Hermitian matrices...

We can for sure sort them (perf hit is likely negligible), but where is this assumed?

KlausC · 2023-02-13T22:38:56Z

stdlib/LinearAlgebra/src/eigen.jl

 ```
 """
 struct Eigen{T,V,S<:AbstractMatrix,U<:AbstractVector} <: Factorization{T}
    values::U
    vectors::S
-    Eigen{T,V,S,U}(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V,S,U} =
-        new(values, vectors)
+    are_vectors_unitary::Bool


This new field will not be required with #48657, which goes a step further in storing the left eigenvectors. For the unitary case F.vectors === F.vectorsl will have th the same information.

stdlib/LinearAlgebra/src/eigen.jl

stdlib/LinearAlgebra/src/symmetriceigen.jl

stdlib/LinearAlgebra/src/tridiag.jl

nanosoldier · 2023-02-13T23:12:02Z

Your package evaluation job has completed - possible new issues were detected.
A full report can be found here.

antoine-levitt · 2023-02-14T06:19:44Z

Doesn't look like there's anything wrong with the pkgeval job

dkarrasch · 2023-02-14T09:44:45Z

but where is this assumed?

I think I confused it with SVD and the default ordering of eigenvalues.

KlausC · 2023-02-14T10:41:29Z

Another issue is that eigenvalues are no longer ordered in ascending order, as we typically assume for values, certainly for symmetric/Hermitian matrices...

We can for sure sort them (perf hit is likely negligible), but where is this assumed?

This is in ?eigen:

By default, the eigenvalues and vectors are sorted lexicographically by (real(λ),imag(λ)). A different comparison function by(λ) can be passed to
  sortby, or you can pass sortby=nothing to leave the eigenvalues in an arbitrary order. Some special matrix types (e.g. Diagonal or SymTridiagonal)
  may implement their own sorting convention and not accept a sortby keyword.

KlausC · 2023-02-14T13:45:52Z

stdlib/LinearAlgebra/src/eigen.jl

+det(A::Eigen) = prod(A.values)
+tr(E::Eigen) = sum(E.values)
+
+*(E::Eigen, x::AbstractVector) = E.vectors * (E.values .* (E.are_vectors_unitary ? (E.vectors'x) : (E.vectors \ x)))


That means a O(n^3) operation for each call. It would make sense to store inv(vectors)' in the Eigen structure here ( and remove are_vectors_unitary). That is also valid for the next function and all others where inv(vectors) or similar are called on the fly.

Yes but this pays the cost upfront, even if the user doesn't need to do anything except computing the values/vectors, right?

Yes but this pays the cost upfront

Cannot understand this argument. If you want to multiply the same E with two different vectors, and you don't know the second one before you have evaluated the first one, you have to E.vectors \ v twice, which means two LU-factorizations of E.vectors. That is what I meant.

Oh, you mean in this function, compute the inverse of the vectors, and mutate E? I don't think we should do something as sophisticated as this, and having Eigen be mutable seems like begging for trouble. These are mostly convenience methods for playing around at the REPL or quick scripts, they are not meant to be used for performance-sensitive work, so over-engineering this adds complexity for not much gain.

KlausC · 2023-02-14T13:57:30Z

stdlib/LinearAlgebra/src/eigen.jl

+mul!(x, E::Eigen, y::AbstractVector) = copyto!(x, E*y)
+
+inv(E::Eigen) = Eigen(map(inv, E.values), copy(E.vectors), E.are_vectors_unitary)
+for f in _matrix_functions


The restriction of applicable functions to this list is a burden. What about adding an additional function
like apply(f, E) = ... to allow for user defined functions?
remark:
unfortuantely (f::Function)(E) = ... is not supported by Julia.

I'm concerned about exporting a function here, esp since the functionality is relatively half-baked (it doesn't apply to general matrices, just to Eigen). I think it would make a lot more sense to have a MatrixFunctions package that does this properly (including for non-hermitian matrices)

it doesn't apply to general matrices, just to Eigen

It would be the typical workflow to first evaluate eigen(A), then f1(A), f2(A), ... . Again it would not make sense to call eigen(A) for each individual evaluation of f(A). That would just be a convenience, if you know in advance, that the eigen-decomposition is only needed for a single function evaluation.
Also your argument is valid for sin, cosh, exp, and brothers in the same way as for apply, so I don't understand your concern.

This operation is valid for all diagonalizable matrices, not only Hermitian or Symmetric, or normal ones.
In the case of multiple eigenvalues and non-diagonalizable matrices, eigen cannot help much.

Also your argument is valid for sin, cosh, exp, and brothers in the same way as for apply, so I don't understand your concern.

Adding methods to these functions does not add API surface

This operation is valid for all diagonalizable matrices, not only Hermitian or Symmetric, or normal ones.

Technically, yes. In practice, using it for strongly non-normal matrices is a bad idea (see eg the classic paper "N bad ways to compute matrix exponentials" or whatever it's called), and there are better ways to do it. Ideally, a MatrixFunctions package would deal with this. This PR is just an opportunity to define methods that ought to exist already.

antoine-levitt · 2023-03-15T15:18:50Z

Can we merge this? CI appears to be OK except for some unrelated problems

dkarrasch · 2023-03-16T11:02:12Z

stdlib/LinearAlgebra/src/eigen.jl

+*(x::Adjoint{<:Any, <:AbstractVector}, E::Eigen) = E.are_vectors_unitary ?
+                                                   ((x*(E.vectors)) * Diagonal(E.values)) * (E.vectors') :
+                                                   ((x*(E.vectors)) * Diagonal(E.values)) / E.vectors
+mul!(x, E::Eigen, y::AbstractVector) = copyto!(x, E*y)


What's the point of providing this method if it doesn't provide any improvement over E*y(and makes it actually even worse, against expectations)?

Very reasonable point, I'll remove it!

dkarrasch

Looking at the matrix functions, I don't think it makes sense to reorder the inverted eigenvalues (and the corresponding eigenvectors). For consistency, we would need to do the same for those Eigen objects, which seems awkward. But then perhaps it would be good to document that after applying inv or matrix functions, you can no longer assume that the first/last entry of the values field contains the smallest/largest eigenvalue by default.

antoine-levitt · 2023-03-22T13:37:48Z

OK, I added a comment on that

antoine-levitt added 3 commits February 6, 2023 09:38

Factor out matrix functions

8ba0f1f

Store the information that the eigenvectors are unitary

6068e32

Extend the Eigen API

b994e3c

antoine-levitt added the linear algebra Linear algebra label Feb 6, 2023

Fix bugs

a53a04b

news + docs

f7d8239

dkarrasch approved these changes Feb 9, 2023

View reviewed changes

garrison reviewed Feb 12, 2023

View reviewed changes

KlausC reviewed Feb 13, 2023

View reviewed changes

stdlib/LinearAlgebra/src/eigen.jl Outdated Show resolved Hide resolved

KlausC reviewed Feb 13, 2023

View reviewed changes

stdlib/LinearAlgebra/src/symmetriceigen.jl Show resolved Hide resolved

KlausC reviewed Feb 13, 2023

View reviewed changes

stdlib/LinearAlgebra/src/tridiag.jl Show resolved Hide resolved

KlausC mentioned this pull request Feb 13, 2023

LinearAlgebra: extend eigen with left eigenvectors and errorbounds #48657

Open

KlausC reviewed Feb 14, 2023

View reviewed changes

Merge branch 'master' into eigen_dispatches

478b8b8

dkarrasch reviewed Mar 16, 2023

View reviewed changes

review

348fd04

dkarrasch approved these changes Mar 20, 2023

View reviewed changes

Add note on ordering

c820151

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Taking Eigen seriously #48549

Taking Eigen seriously #48549

antoine-levitt commented Feb 6, 2023 •

edited

dkarrasch commented Feb 6, 2023

antoine-levitt commented Feb 6, 2023

dkarrasch commented Feb 6, 2023

antoine-levitt commented Feb 7, 2023

RalphAS commented Feb 8, 2023

antoine-levitt commented Feb 9, 2023

dkarrasch left a comment

antoine-levitt commented Feb 9, 2023

garrison Feb 12, 2023 •

edited

dkarrasch Feb 12, 2023

antoine-levitt Feb 13, 2023

LilithHafner Feb 13, 2023

dkarrasch commented Feb 12, 2023

antoine-levitt commented Feb 13, 2023

KlausC Feb 13, 2023

nanosoldier commented Feb 13, 2023

antoine-levitt commented Feb 14, 2023

dkarrasch commented Feb 14, 2023

KlausC commented Feb 14, 2023

KlausC Feb 14, 2023 •

edited

antoine-levitt Mar 3, 2023

KlausC Mar 3, 2023

antoine-levitt Mar 3, 2023

KlausC Feb 14, 2023

antoine-levitt Mar 3, 2023

KlausC Mar 3, 2023

antoine-levitt Mar 3, 2023

antoine-levitt commented Mar 15, 2023

dkarrasch Mar 16, 2023

antoine-levitt Mar 16, 2023

dkarrasch left a comment

antoine-levitt commented Mar 22, 2023

Taking Eigen seriously #48549

Are you sure you want to change the base?

Taking Eigen seriously #48549

Conversation

antoine-levitt commented Feb 6, 2023 • edited

dkarrasch commented Feb 6, 2023

antoine-levitt commented Feb 6, 2023

dkarrasch commented Feb 6, 2023

antoine-levitt commented Feb 7, 2023

RalphAS commented Feb 8, 2023

antoine-levitt commented Feb 9, 2023

dkarrasch left a comment

Choose a reason for hiding this comment

antoine-levitt commented Feb 9, 2023

garrison Feb 12, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkarrasch commented Feb 12, 2023

antoine-levitt commented Feb 13, 2023

Choose a reason for hiding this comment

nanosoldier commented Feb 13, 2023

antoine-levitt commented Feb 14, 2023

dkarrasch commented Feb 14, 2023

KlausC commented Feb 14, 2023

KlausC Feb 14, 2023 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

antoine-levitt commented Mar 15, 2023

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkarrasch left a comment

Choose a reason for hiding this comment

antoine-levitt commented Mar 22, 2023

antoine-levitt commented Feb 6, 2023 •

edited

garrison Feb 12, 2023 •

edited

KlausC Feb 14, 2023 •

edited