Julia: Justification for eig, eigvecs, eigvals and eigfact?

[nalimilan]

I wonder whether we really need to offer eig, eigvecs and eigfact at the same time. eig and eigvecs are respectively a two- and a one-line wrapper around eigfact.

Also, maybe eigvals! could be replaced with an eigfact! function taking the vector in which to store eigenvectors. This would allow getting rid of eigvals, which is yet another wrapper around eigfact.

The same applies to the svd* family of functions.

A more speculative note: if tuple destructuring was extended to immutables, one could write directly val, vec = eigfact(X). This would offer a perfect replacement for eig.

[tkelman]

I think it could be worth cleaning up some namespace here.

[andreasnoack]

The reason to define eigvecs is to support inverse iteration when the values are known. Maybe this could be merged into a method for eigfact.

As I read your proposal, eigfact wouldn't return a Factorization anymore which I think would be unfortunate.

Right now, a problem with the factorizations is that extracting the factors is not type stable. Hopefully, something like dot overloading will fix that and cleaning up the namespace will probably be easier at that point.

[nalimilan]

The reason to define eigvecs is to support inverse iteration when the values are known. Maybe this could be merged into a method for eigfact.

Indeed. Maybe with a keyword argument?

As I read your proposal, eigfact wouldn't return a Factorization anymore which I think would be unfortunate.

No, my proposal is to keep only the variant returning a Factorization object.

Right now, a problem with the factorizations is that extracting the factors is not type stable. Hopefully, something like dot overloading will fix that and cleaning up the namespace will probably be easier at that point.

Maybe @pure could help with this: https://github.com/JuliaStats/DataFrames.jl/issues/744#issuecomment-195854129 I've played a bit with it, and it doesn't seem to work here, but maybe it could be improved.

[nalimilan]

@vtjnash Do you think @pure could allow making getindex(::Eigen, ::Symbol) type-stable when passed a literal symbol (see comment above)?

[vtjnash]

if tuple destructuring was extended to immutables

what makes you think that it isn't? all it needs is an iterable:

julia> a,b = (1=>2)
1=>2

julia> a
1

julia> b
2

Do you think @pure could allow making getindex(::Eigen, ::Symbol) type-stable when passed a literal symbol (see comment above)?

that's an interesting combination. i don't see an obvious way to make this work however. i like the iteration version. i think carnaval has been talking about providing a way to make the getindex version easier as well, perhaps in v0.6.

[nalimilan]

@vtjnash I guess I should have tried destructuring a non-tuple. That's cool.

Regarding @pure, how is it that it works for https://github.com/JuliaStats/DataFrames.jl/issues/744#issuecomment-195854129, but not here? Anyway, we can wait for 0.6.

[jiahao]

eigvals, which is yet another wrapper around eigfact.

This statement is correct only as the generic fallback for types like ordinary matrices. For some specialized matrix types (notably SymTridiagonal) eigvals can be computed with less effort than the full eigfact.

I introduced eigvecs as a workaround for the inherent type instability in eigfact(A)[:vectors] and I got tired of typing eigfact(A)[:vectors]::TheCorrectMatrixTypeOfEigenvectors constantly in my code. Until we fix the eigfact API I really don't think it's a good idea to remove eigvecs.

eig is really a Matlab compatibility function and I'd be fine with getting rid of it. I suspect that many users would not, however.

[nalimilan]

eig could go to MATLAB.jl, like other compatibility functions.

But it doesn't look like we're ready to replace it with eigfact. Indeed, I was right that destructuring an Eigen object doesn't work: it works for any iterable, but Eigen isn't iterable, and it doesn't look like it should be. It would be nice to be able to destructure any object, but this conflicts with the definition that destructuring uses the iteration protocol.

[vtjnash]

Indeed, I was right that destructuring an Eigen object doesn't work

i was roughly thinking that the Eigen return type would be replaced by a (LazyEigval, LazyEigvec). but i'm sure the details of such a proposal would need to be worked out further first.

[nalimilan]

i was roughly thinking that the Eigen return type would be replaced by a (LazyEigval, LazyEigvec). but i'm sure the details of such a proposal would need to be worked out further first.

Hmm, I'm not sure how that would work. Returning an immutable sounds quite natural here (if it weren't for the impossibility of destructuring it).

I've played a bit more with @pure and Val, with no luck. Though it seems the problem isn't with either of these, as Val{:values} is correctly inferred: it seems that && branches are not removed by type inference (which I guess is a known issue). Example:

function getindex2{T}(x, ::Type{T})
    T <: Val{:values} && return x.values
    T <: Val{:vectors} && return x.vectors
    throw(KeyError(T))
end

function f(x::Symbol)
    Base.@_pure_meta
    Val{x}
end

function g()
    A = Symmetric([1 0; 0 1])
    getindex2(eigfact(A), f(:values))
end

@code_warntype g()

UPDATE: getindex2 doesn't need to be @pure.
UPDATE2: using a Symmetric array is needed to get an inferrable result even with eig

[vtjnash]

Since getindex2 is not particularly const, it's not a good idea to mark it. This transformation is also not what @pure is supposed to be for (which is to allow memoization of the results of constant expressions).

[nalimilan]

Since getindex2 is not particularly const, it's not a good idea to mark it. This transformation is also not what @pure is supposed to be for (which is to allow memoization of the results of constant expressions).

Actually, I don't think getindex2 needs to be declared as @pure, right? (At least, in my tests it doesn't change anything.) This was just an extreme attempt at trying to get inference to work. The key here is f().

Is it expected that the T <: Val{:values} && branch isn't eliminated at compile time?

[carnaval]

Is it expected that the T <: Val{:values} && branch isn't eliminated at compile time?

it should be eliminated unfortunately this happens at codegen time and it's too late to get the inferred type right. This is indeed a known issue.

[nalimilan]

OK, thanks. So it looks like we'll have to wait for this to work.

[nalimilan]

Looks like inference has progressed in 0.5 and master (compared to 0.4.7). Now, the example from my comment above is correctly inferred. Unfortunately, I've not been able to make the destructuring using iteration protocol type stable. In the two following functions, the type of the first variable (a) is correctly inferrred, but not of the second one (b). Sounds promising, though; maybe there's just a small piece missing to get this to work.

import Base.Eigen
Base.start(::Eigen) = Val{:values}
Base.done(::Eigen, state) = state <: Val{:done}
function Base.next{T}(x::Eigen, state::Type{T})
    if T <: Val{:values}
        return x.values, Val{:vectors}
    elseif T <: Val{:vectors}
        return x.vectors, Val{:done}
    else
        throw(KeyError(T))
    end
end

function f(A)
    a, b = eigfact(A)
    a, b
end

A = Symmetric([1 0; 0 1])
@code_warntype f(A)

function g(A)
    f = eigfact(A)
    s = start(f)
    a, s = next(f, s)
    b, s = next(f, s)
    a, b
end

@code_warntype g(A)

[nalimilan]

For triage: is there anything we could be done for 0.7 here? If we expect to make eigfact(A)[:vectors] type-stable (maybe thanks to IPO?), it would make sense to clean up the API now and deprecate eig before it's too late.

[ViralBShah]

If we move eig to Matlab compatibility, we should do it for all dense linalg.

[JeffBezanson]

Triage thinks we should wait until a solution for eigfact(A)[:vectors] exists, and leave this alone until then.

[nalimilan]

@andreasnoack Do you plan to change any of the functions mentioned here?

[andreasnoack]

Yes. I'll try to continue the work in #25187 and see if can work generally.

[Sacha0]

In a similar vein, is there any reason for chol to exist alongside cholfact? The few methods for the former seem like they could be subsumed into the latter? Best!

[ViralBShah]

+1 to cleaning up things that set us up for the long term.

[andreasnoack]

Fixed by https://github.com/JuliaLang/julia/pull/27212