Address type stability issues in #12574 and fix a bug or two #12594

kshyatt · 2015-08-13T01:03:45Z

The trix methods should have been trix! since they return references. The methods for Unit{Upper/Lower]Triangular matrices were wrong (now fixed) and I also caught a bug in istriu/istril for Tridiagonals.

cc: @andreasnoack and @tkelman

tkelman · 2015-08-13T01:19:53Z

base/linalg/bidiag.jl

@@ -108,7 +108,7 @@ ctranspose(M::Bidiagonal) = Bidiagonal(conj(M.dv), conj(M.ev), !M.isupper)
 istriu(M::Bidiagonal) = M.isupper || all(M.ev .== 0)
 istril(M::Bidiagonal) = !M.isupper || all(M.ev .== 0)

-function tril(M::Bidiagonal, k::Integer=0)
+function tril!(M::Bidiagonal, k::Integer=0)


this could use fill! instead of allocating new zeros arrays, right?

I tested this locally and it works. Will wait for CI to finish then update.

I think this applies throughout for (most of?) the rest of the types too. Anywhere it's possible to make the ! versions allocation-free and in-place, then I think that would be a good goal. And aim for type-stability but with the smallest amount of type widening that makes sense. The one-arg case could possibly be made to return a more specialized type than the two-arg case for some types since it's more specific in its behavior, but that would require extra methods rather than using default argument values and could be left as a future enhancement.

kshyatt · 2015-08-13T03:30:29Z

Updated! Should be far fewer allocs now.

tkelman · 2015-08-13T03:38:05Z

base/linalg/diagonal.jl

+    elseif k == 0
+        return D
+    else
+        return zeros(D)


fill!(D.diag, 0) here

Noooo ok and I need a whitespace fix. This is not my night.

Not a big deal, stop me if I'm getting too nitpicky.

No it's good. I want this to work and work well.

tkelman · 2015-08-13T03:56:50Z

Not this PR's fault, but are the docs for tril(M, k) wrong? It should be k-th superdiagonal, same as triu(M, k), not k-th subdiagonal like the docs currently say, right?

kshyatt · 2015-08-13T03:57:55Z

@jiahao wanted that phrasing, ask him :)

tkelman · 2015-08-13T04:00:53Z

base/linalg/triangular.jl

    elseif k == 0
        return UnitUpperTriangular(eye(A))
    else
-        return UnitUpperTriangular(triu(tril(A.data,k)))
+        return UnitUpperTriangular(triu!(tril!(A.data,k)))


this one hurts my head a little - maybe the easiest way to get this right, and get the rest of tril!(A::UnitUpperTriangular, k) for almost free would be first set the diagonals of A.data to one, then return tril!(UpperTriangular(A.data), k) ?

tkelman · 2015-08-13T04:05:06Z

But that would be the opposite meaning. What the docs describe is not the implemented behavior:

julia> A = rand(3,3)
3x3 Array{Float64,2}:
 0.560414  0.632475  0.901667
 0.342709  0.46708   0.133265
 0.604388  0.354852  0.706773

julia> tril(A, 1)
3x3 Array{Float64,2}:
 0.560414  0.632475  0.0
 0.342709  0.46708   0.133265
 0.604388  0.354852  0.706773

julia> triu(A, 1)
3x3 Array{Float64,2}:
 0.0  0.632475  0.901667
 0.0  0.0       0.133265
 0.0  0.0       0.0

tkelman · 2015-08-13T04:09:56Z

base/linalg/triangular.jl

-tril(A::UnitLowerTriangular,k::Integer=0) = UnitLowerTriangular(tril(tril(A.data),k))
+tril!(A::UnitLowerTriangular,k::Integer=0) = UnitLowerTriangular(tril!(tril!(A.data),k))
+
+tril(A::UpperTriangular,k::Integer=0)     = tril!(copy(A),k)


should the generic fallbacks here

julia/base/linalg/generic.jl

Lines 41 to 44 in 24a92a9

triu(M::AbstractMatrix) = triu!(copy(M))

tril(M::AbstractMatrix) = tril!(copy(M))

triu!(M::AbstractMatrix) = triu!(M,0)

tril!(M::AbstractMatrix) = tril!(M,0)

be tweaked to cover all of these?

tkelman · 2015-08-13T04:42:41Z

This is turning out to be a surprisingly complicated set of methods to get right. Someone who's good at diagrams could make a pretty neat picture about the algebraic closedness of the different operations here for all types and super/sub diagonal cases.

kshyatt · 2015-08-13T04:44:52Z

OK I think I have something semi-reasonable now. Let's give it a whirl.

tkelman · 2015-08-13T04:48:26Z

base/linalg/diagonal.jl

+    n = size(D,1)
+    if abs(k) > n
+        throw(ArgumentError("requested diagonal, $k, out of bounds in matrix of size ($n,$n)"))
+    elseif k != 0


only if > 0, right?

~~Can't k be negative?~~ Oh, derp, I see.

jiahao · 2015-08-13T04:52:10Z

Someone who's good at diagrams could make a pretty neat picture about the algebraic closedness of the different operations here for all types and super/sub diagonal cases.

See #8240 where I discuss how the algebraic structure shows up as banded matrices.

It should be k-th superdiagonal, same as triu(M, k), not k-th subdiagonal like the docs currently say, right?

I don't recall... but certainly tril(randn(4,4), -1) produces the matrix below and starting from the 1st subdiagonal,

tkelman · 2015-08-13T04:55:41Z

base/linalg/triangular.jl

+        for i in diagind(A)
+            A.data[i] = one(eltype(A))
+        end
+        return UpperTriangular(triu!(tril!(A.data,k)))


triu! still not necessary here

tkelman · 2015-08-13T04:58:17Z

tril(randn(4,4), -1) produces the matrix below and starting from the 1st subdiagonal,

Right, so either the docs should say "-k th" subdiagonal, or "k th superdiagonal."

jiahao · 2015-08-13T05:03:07Z

yes

tkelman · 2015-08-13T05:14:17Z

the algebraic structure shows up as banded matrices.

Well, banded has a much nicer, easier to deal with algebraic structure than our current menagerie. The less uniform more complicated version we have now would make for a more interesting picture, but that's not necessarily an endorsement.

tkelman · 2015-08-13T05:30:57Z

base/linalg/triangular.jl

 end

-tril(A::UnitLowerTriangular,k::Integer=0) = UnitLowerTriangular(tril(tril(A.data),k))
+tril(A::UpperTriangular,k::Integer=0)     = tril!(copy(A),k)


Not sure why github collapsed #12594 (comment), but that still applies. I think the others that have been collapsed as of df4c774 have been addressed, but will review again tomorrow.

edit: links to collapsed comments apparently don't always work so well, that was

should the generic fallbacks here

julia/base/linalg/generic.jl

Lines 41 to 44 in 24a92a9

triu(M::AbstractMatrix) = triu!(copy(M))

tril(M::AbstractMatrix) = tril!(copy(M))

triu!(M::AbstractMatrix) = triu!(M,0)

tril!(M::AbstractMatrix) = tril!(M,0)

be tweaked to cover all of these?

stevengj · 2015-08-13T12:22:17Z

base/linalg/symmetric.jl

-triu(A::Symmetric,k::Integer=0) = triu(A.data,k)
+function tril(A::Hermitian, k::Integer=0)
+    if A.uplo == 'U' && k <= 0
+        return tril(A.data',k)


Shouldn't this be tril!(A.data',k), since the data' already makes a copy? Similarly below.

Good catch. I've updated the PR to fix this.

StefanKarpinski · 2015-08-13T13:39:13Z

This is turning out to be a surprisingly complicated set of methods to get right.

That's an excellent sign of something worth having in a library :-)

kshyatt · 2015-08-13T20:35:38Z

I have updated the PR to fix @tkelman and @stevengj's catches. Does this look ok, now?

tkelman · 2015-08-13T20:42:23Z

base/linalg/symmetric.jl

-triu(A::Hermitian,k::Integer=0) = triu(A.data,k)
-tril(A::Symmetric,k::Integer=0) = tril(A.data,k)
-triu(A::Symmetric,k::Integer=0) = triu(A.data,k)
+function tril(A::Hermitian, k::Integer=0)


This could be Union{Hermitian,Symmetric} right? The implementations look like exact copies.

It passes tests locally with this change. Do we need to run CI again?

tkelman · 2015-08-13T20:49:43Z

Getting there, but not quite. Not yet addressed:

And will want to double-check that comments left so far are fixed for both mirrored sets of methods.

kshyatt · 2015-08-13T21:24:38Z

Ok, I've gone through on my end and dealt with everything but the last one. #12594 (comment). I'm not 100% how to tweak the fallbacks - advice?

tkelman · 2015-08-13T21:30:09Z

The decision to make there is whether adding general fallbacks that look like

triu(M::AbstractMatrix, k::Integer) = triu!(copy(M),k)
tril(M::AbstractMatrix, k::Integer) = tril!(copy(M),k)

makes sense, now that you've implemented a whole bunch of these methods for many of the subtypes of AbstractMatrix. Would be cleaner and should subsume the dozen-ish versions you've added here.

edit: sorry, just the first 2, the 3rd and 4th lines there were meaningless

kshyatt · 2015-08-13T21:40:10Z

Shall I give a go and report back?

E: It's working!!!

tkelman · 2015-08-13T22:34:58Z

base/linalg/symmetric.jl

-triu(A::Hermitian,k::Integer=0) = triu(A.data,k)
-tril(A::Symmetric,k::Integer=0) = tril(A.data,k)
-triu(A::Symmetric,k::Integer=0) = triu(A.data,k)
+function tril(A::Union{Hermitian,Symmetric}, k::Integer=0)


Ack, I was wrong, while the earlier implementations were identical for Hermitian vs Symmetric, they actually shouldn't be. We can have Symmetric for a complex element type in which case these should be using .' instead of '. But for Hermitian they should use '.

oh nooooo!

Ok, let CI run, then update?

Address type stability issues in #12574 and fix a bug or two

ScottPJones · 2015-08-14T17:30:57Z

I noticed when reviewing the previous change for error messages and coverage, and looking at this one, that the BunchKaufman, Cholesky, CholeskyPivoted, Symmetric, Hermitian types use uplo::Char, but the Bidiagonal type uses a isupper::Bool.
I think the code in general would be more efficient if it consistently used isupper::Bool, so I wonder if I'm missing something.

tkelman · 2015-08-14T22:25:19Z

That inconsistency is a bit annoying, but they serve slightly different purposes. For the ~~factorization and~~ Symmetric/Hermitian objects, the matrices aren't really completely upper or lower triangular - the uplo parameter indicates which triangle the important information is stored in. Using uplo::Char also maps more directly to the BLAS and LAPACK API's. I'm not sure if we are binding to much of LAPACK for Bidiagonal, but it might make more sense to switch Bidiagonal from using isupper::Bool to uplo::Char for consistency.

ScottPJones · 2015-08-15T13:07:04Z

I had looked at the bindings, but even so, I think overall it would be more efficient if all of them were changed to use isupper::Bool (I can attempt a PR later). For the bindings, it's only a single instruction (on x86 at least) to convert the bool to 'L'/'U'. i.e. leal 76(%edi,%edi,8), %eax, and it can save 3 bytes in the structures 😀, besides simplifying a lot of places in the code I looked at.

tkelman · 2015-08-16T02:22:46Z

Remember these are Fortran API's, so you have to pass a reference to the char anyway. I doubt such a PR would be especially popular with the people who've been developing and maintaining the blas and lapack bindings.

ScottPJones · 2015-08-16T02:28:41Z

I wasn't aware of that for the Fortran interface. Why do you believe something that could simplify the code would not be popular with them?
It would still be just an instruction to pass either a reference to a L or U.

tkelman · 2015-08-16T02:30:50Z

I'm guessing, but it seems like unnecessary churn that would make our code map less directly to the standard-for-decades, widely used API's for linear algebra. You can always open a PR, but "simplify the code" is a relative term that is going to mean different things to different people.

ScottPJones · 2015-08-16T05:41:57Z

There may need to be some "churn", it looks like the code as is won't even work on a big-endian machine. It's passing a pointer to a 32-bit Char value to something expecting a pointer to a 8-bit char.
It will get '\0' on big-endian systems, not 'L' or 'U'.

stevengj · 2015-08-17T15:30:04Z

@ScottPJones, if you pass &uplo::Char in ccall for a Ptr{UInt8} argument, Julia will do the right thing: it converts the Char to UInt8 before passing the address. And this conversion is always safe because the Char is checked to have an ASCII value before the conversion.

tkelman · 2015-08-17T20:52:23Z

x-ref https://groups.google.com/forum/#!topic/julia-dev/MMsJvqIJ3yY

kshyatt added the linear algebra Linear algebra label Aug 13, 2015

tkelman reviewed Aug 13, 2015
View reviewed changes

kshyatt force-pushed the ksh/fixtrix branch from 9a566cb to 4514ac8 Compare August 13, 2015 03:30

tkelman reviewed Aug 13, 2015
View reviewed changes

kshyatt force-pushed the ksh/fixtrix branch from 4514ac8 to df4c774 Compare August 13, 2015 04:45

tkelman reviewed Aug 13, 2015
View reviewed changes

kshyatt force-pushed the ksh/fixtrix branch from df4c774 to 161938c Compare August 13, 2015 07:19

stevengj reviewed Aug 13, 2015
View reviewed changes

kshyatt added 4 commits August 13, 2015 11:04

fixed up trix for triangular

dd009d6

fixed tridiag

7465943

fixed diagonal and added istrix

4d50e23

Fixed up bidiag too

095e0c4

kshyatt force-pushed the ksh/fixtrix branch from 161938c to a4f616a Compare August 13, 2015 18:04

tkelman reviewed Aug 13, 2015
View reviewed changes

kshyatt force-pushed the ksh/fixtrix branch from a4f616a to a67c4e2 Compare August 13, 2015 22:14

tkelman reviewed Aug 13, 2015
View reviewed changes

kshyatt added 2 commits August 13, 2015 18:55

Fixed symmetric

1c08f8d

Fix residual issues and add fallback methods

851e7a1

kshyatt force-pushed the ksh/fixtrix branch from a67c4e2 to 851e7a1 Compare August 14, 2015 02:37

kshyatt added a commit that referenced this pull request Aug 14, 2015

Merge pull request #12594 from JuliaLang/ksh/fixtrix

b59b0dd

Address type stability issues in #12574 and fix a bug or two

kshyatt merged commit b59b0dd into master Aug 14, 2015

kshyatt deleted the ksh/fixtrix branch August 14, 2015 08:12

	triu(M::AbstractMatrix) = triu!(copy(M))
	tril(M::AbstractMatrix) = tril!(copy(M))
	triu!(M::AbstractMatrix) = triu!(M,0)
	tril!(M::AbstractMatrix) = tril!(M,0)

Address type stability issues in #12574 and fix a bug or two #12594

Address type stability issues in #12574 and fix a bug or two #12594

Conversation

kshyatt commented Aug 13, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

kshyatt commented Aug 13, 2015

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

tkelman commented Aug 13, 2015

kshyatt commented Aug 13, 2015

Choose a reason for hiding this comment

tkelman commented Aug 13, 2015

Choose a reason for hiding this comment

tkelman commented Aug 13, 2015

kshyatt commented Aug 13, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

jiahao commented Aug 13, 2015

Choose a reason for hiding this comment

tkelman commented Aug 13, 2015

jiahao commented Aug 13, 2015

tkelman commented Aug 13, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

StefanKarpinski commented Aug 13, 2015

kshyatt commented Aug 13, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkelman commented Aug 13, 2015

kshyatt commented Aug 13, 2015

tkelman commented Aug 13, 2015

kshyatt commented Aug 13, 2015

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ScottPJones commented Aug 14, 2015

tkelman commented Aug 14, 2015

ScottPJones commented Aug 15, 2015

tkelman commented Aug 16, 2015

ScottPJones commented Aug 16, 2015

tkelman commented Aug 16, 2015

ScottPJones commented Aug 16, 2015

stevengj commented Aug 17, 2015

tkelman commented Aug 17, 2015