type-instability in LinAlg factorizations getindex declarations

[vtjnash]

factorizations.jl has a number of getindex(x, ::Symbol) functions, all of whose return types are dependent upon the value of the symbol. This seems to be a performance trap, no? (for example: julia can't infer the return type of eig)

[ViralBShah]

I believe we had discussed this in an earlier issue, and @JeffBezanson had said this should be fine.

[JeffBezanson]

If it is considered a bug every time an inferred type isn't concrete, then
julia should be statically typed.
On May 5, 2013 10:47 AM, "Viral B. Shah" [email protected] wrote:

I believe we had discussed this in an earlier issue, and @JeffBezansonhttps://github.com/JeffBezansonhad said this should be fine.

—
Reply to this email directly or view it on GitHubhttps://github.com//issues/12
.

[JeffBezanson]

JuliaLang/julia#1974 will fix this.

[ViralBShah]

This hasn't been an issue for a long time.

[andreasnoack]

I'd actually like fix this issue as it can cause significant slow down as we saw recently for the correlation benchmark. The solution is probably to use typed instead of symbols for the indexing and declaration of special matrix and factorizations.

[ViralBShah]

Ok, I missed that. Is there a relevant issue we can link here?

[andreasnoack]

We discussed the Triangular case in #144

[ViralBShah]

Thanks!

[andreasnoack]

@jiahao is giving some more context and examples in #241. A couple of solutions have been mentioned

1. The Val solution, e.g. qrfact(A)[Val{:Q}]
2. Dummy types, e.g. Q, R, P etc. types
3. Accessor functions like getQ which is already defined, but unexported in LinAlg.
4. Same as 3. but implemented as . overloaded methods.

I think 1. is simply too ugly and I believe/hope that eventually we'll be able to fix this in a better way. 2. reseves too many one letter names. 3. defines way too many functions each with a single simple use case. 4. might be the best bet but depends on JuliaLang/julia#1974. Value specialization might also help us. @carnaval can tell if that is true or if I'm missing the point.

[vtjnash]

defines way too many functions each with a single simple use case.

i think you are underselling the value of having one function to do one thing.

[andreasnoack]

It would be something like getQ, getR, getP, getp, getValues, getVectors, getLeftVectors, getRightVectors, getT, getL, getU, getC, getD, getU, getV.

[jiahao]

JuliaLang/julia#12747 introduced a method of svdvals that serves as an accessor into the singular values in an SVD object. Similar methods for eigvals and eigvecs can be defined.

The motivating example in #241 is stripped down from IterativeSolvers.svdvals_gkl in JuliaLinearAlgebra/IterativeSolvers.jl@a7484fa.

The code basically does this

function svdvals_gkl(A, threshold)
   T=typeof(real(one(eltype(A)))
   otherstate = initialize(A)
   converged_values = T[] #<--
   while !enough(converged_values)
    M, otherstate = make_small_matrix(A, otherstate, converged_values) #<--
    S = svdfact(M)
    svals = S[:S] #<--
    svalerrs = compute_errors(S)
    converged_values = svals[svalerrs .< threshold] #<--
  end
end

In the full code, the fact that make_small_matrix_from depends on data derived from svals confuses type inference enough that at some point converged_values is inferred to be of type Any, which ruins performance.

[carnaval]

The long term elegant solution is probably the getfield overloading (which is the same as specializing on symbol values).
In the meantime maybe split it into several functions depending roughly on return types ?
matrix_factor(fact,:L) (better name needed).

[carnaval]

Actually you can probably do it this way : F[Matrix,:L]
where using the wrong type for a symbol give a nice error ?

There may even be some symbols that you would want to return different things depending on the type ?

[jiahao]

That could be one solution.

There may even be some symbols that you would want to return different things depending on the type ?

Possibly. You might want F[Vector, :values] to return eigenvalues as a list but F[Diagonal, :values] to return them as a diagonal matrix.

Less trivial is QR[Matrix, :Q] == full(QR[:Q]) as opposed to the fancy QR objects that can come out of qrfact.

[StefanKarpinski]

Given the lack of progress at this point, this is not going to make 0.6.

[StefanKarpinski]

Bump?

[andreasnoack]

We are waiting for a decision on the dot overload issue. Accessing the factors could be made convenient and type stable if we had some kind of dot overloading. E.g. something like F.L() instead of F[:L]. If we don't get dot overloading or some kind of type local functions/types, we'd probably have to go the Val even though I'd like to avoid it. Defining normal generic get functions as described in #12 doesn't seem desirable.

[vtjnash]

Close? The original issues I brought up aren't valid. It seems like F[:L] is the best API right now. When we add getfield-overloading, we can consider updating it to add F.L as syntax for that.

[StefanKarpinski]

The original issue you brought up still seems valid to me. We need either enough constant propagation to infer the type of F[:L] based on the constant symbol argument or we need dot overloading, neither of which we have. I agree, however, that this is not 1.0-blocking since we can add either of those in a 1.x release.

[andreasnoack]

Yes. Let's keep this open until we can extract factors in a type stable way.

[KristofferC]

Constant propagation + getfield overloading seems to close this 🎉. Reopen if not.