Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added ConjArray wrapper type for conjugate array views #20047

Merged
merged 1 commit into from
Feb 9, 2017

Conversation

andyferris
Copy link
Member

By default, this is used only in conjugation with RowVector, so that both transpose(vec) and ctranspose(vec) both return views.

However, this will pave the way for views for transposed matrices and replacement of Ac_mul_B, etc by specialized methods on *, simplifying parsing.

@andyferris
Copy link
Member Author

Follow-up of #19670

@propagate_inbounds getindex{T,N}(a::ConjArray{T,N}, i::Int) = conj(getindex(a.parent, i))
@propagate_inbounds getindex{T,N}(a::ConjArray{T,N}, i::Vararg{Int,N}) = conj(getindex(a.parent, i...))
@propagate_inbounds setindex!{T,N}(a::ConjArray{T,N}, v, i::Int) = setindex!(a.parent, conj(v), i)
@propagate_inbounds setindex!{T,N}(a::ConjArray{T,N}, v, i::Vararg{Int,N}) = setindex!(a.parent, conj(v), i...)
Copy link
Member

@stevengj stevengj Jan 15, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Seems like it should pass through a stride(A,k) and pointer (and unsafe_convert to a pointer) call to the underlying array. (Similarly for RowVector.)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually, not pointer since the underlying data is not conjugated. But we should provide pointer for RowVector. And we should still provide stride for ConjArray (this is potentially useful for telling you what direction to traverse an array).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, I was hoping to address all the StridedArray behavior at some point...

(it is a bit annoying that it is a type not a trait)

parent::A
end

@inline ConjArray{T,N}(a::AbstractArray{T,N}) = ConjArray{conj_type(T), N, typeof(a)}(a)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We should probably have a similar function that calls similar(parent).

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Agreed

@inline _conj(a::AbstractArray) = ConjArray(a)
@inline _conj{T<:Real}(a::AbstractArray{T}) = a
@inline _conj(a::ConjArray) = parent(a)
@inline _conj{T<:Real}(a::ConjArray{T}) = parent(a)
Copy link
Member

@stevengj stevengj Jan 15, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't think this method is needed?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In functionality, you are right, but this is a ambiguity resolution of the two methods above (not that ConjArray{T} where T<:Real is common, but it is possible for a user to create one, unless you prefer for the constructor to throw in that case?)

@stevengj
Copy link
Member

Should't there be some more specialized linear-algebra operations? e.g. x' * y should turn into a dot call. _conj(a) + _conj(b) should be the same as _conj(a + b) (similarly for *).

@stevengj stevengj added the linear algebra Linear algebra label Jan 15, 2017
NEWS.md Outdated
@@ -185,6 +185,10 @@ Library improvements

* `notify` now returns a count of tasks woken up ([#19841]).

* Introduced a wrapper type for lazy complex conjugation of arrays, `ConjArray`.
Currently, it is used by default for the new `RowVector` type only, and
enforces that both `transpose(vec)` and `ctranspose(vec)` are views not copies.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

copies. -> copies ([#20047]).

# @test which(\, (SparseMatrixCSC, AbstractVecOrMat)).module == Base.SparseArrays
#end
@testset "issue #16548" begin
@test_broken which(\, (SparseMatrixCSC, AbstractVecOrMat)).module == Base.SparseArrays
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this should be replaced with a test that serves the same purpose, if possible, cc @Sacha0

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I agree. I haven't been able to think of an elegant way of doing the same check.

Does anyone know the status/purpose of the "\ revisions"?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Status of \ revisions: Chances are no one has worked on improving the mechanism further since #16548 (particularly given that improving that mechanism depends on compiler improvements, see #16979 (comment)).

Had a brief look at fixing this test. A simple fix would be to retrieve the method matches for \ with (SparseMatrixCSC, AbstractVecOrMat) arguments and check that all such methods live in the SparseArrays module. For example, something along the lines of all(meth -> meth.module == Base.SparseArrays, methods(\, (SparseMatrixCSC, AbstractVecOrMat)).ms) could do the trick. Thoughts? Best!

Copy link
Member Author

@andyferris andyferris Jan 21, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Hmm, OK. At first I thought the issue might have been it going through a rowvector.jl method first before going to the Sparse module, but in fact this may be a disambiguity method that I added to Sparse, so that should work.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

a disambiguity method that I added to Sparse

Correct. That there are now two matching methods rather than one is what caused the existing test to fail:

julia> which(\, (SparseMatrixCSC, AbstractVecOrMat))
ERROR: method match is ambiguous for the specified argument types
Stacktrace:
 [1] which(::Any, ::Any) at ./reflection.jl:715

julia> methods(\, (SparseMatrixCSC, AbstractVecOrMat))
# 2 methods for generic function "\":
\(::SparseMatrixCSC, ::RowVector) in Base.SparseArrays at sparse/linalg.jl:875
\(A::SparseMatrixCSC, B::Union{AbstractArray{T,1},AbstractArray{T,2}} where T) in Base.SparseArrays at sparse/linalg.jl:855

Best!

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I fixed this for now, but I'm guessing that this will need to change again with TransposedMatrix (the thing it is testing may not even be true after that...).

@andyferris
Copy link
Member Author

andyferris commented Jan 20, 2017

Should't there be some more specialized linear-algebra operations? e.g. x' * y should turn into a dot call.

I think this is done now.

_conj(a) + _conj(b) should be the same as _conj(a + b) (similarly for *).

Hmm... I thought that getting rid of the ConjArray wrapper when creating a new array would be desirable. Are you suggesting it's more advantageous to take advantage of BLAS-level 1 routines? Shouldn't Julia generic code be pretty good at adding two vectors? (i.e. I assumed a::AbstractArray + b::AbstractArray would be calculated as map(+, a, b))

EDIT: I believe * should work out fine already due to the changes in #19670.

@andyferris
Copy link
Member Author

andyferris commented Jan 20, 2017

A question: what should (v').' and (v.')' become? Should it give a ConjVector or make a new conjugated copy of v?

The goal here was to make both .' and ' return a view for vectors, but I hadn't defined conj(v::AbstractVector) = ConjArray(v), so this seems a little at odds.

EDIT: actually making a copy would mean it's trivial that the BLAS routines get used in *.

@timholy
Copy link
Member

timholy commented Jan 20, 2017

IIRC, the last time I looked Julia was as fast as BLAS for 1d operations (and indeed faster for small objects because of inlining, which is not possible with BLAS). There may have been some regression in our SIMD capabilities, so that may not be so true anymore, but I'd rather see this implemented in Julia code and then work on our SIMD.

@@ -0,0 +1,23 @@
# This file is a part of Julia. License is MIT: http://julialang.org/license

@testset "RowVector" begin
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

don't think this outer testset is necessary, the test system inserts a test set around each file anyway

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right, thanks for the heads-up. The behavior of Base tests is quite different to what happens with packages, but I had seen an issue (from memory) about using @testset more so I used my usual style and hoped things would be picked up.

(is the hope to unify the package and base approaches in the future?)

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

yes, #18740 and #19567

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Very cool!

@andyferris
Copy link
Member Author

Just saw this. Should we add a milestone for v0.6 for this PR? IMO, it is a (necessary?) refinement of #19670 (but it also adds a feature, so isn't just a bugfix).

@stevengj
Copy link
Member

There are the BLAS3 cases of conjmatrix * conjmatrix, conjtransposematrix * matrix, etcetera, where will it certainly be faster to rewrite in terms of matrix multiplies on the parent matrices.

@stevengj
Copy link
Member

stevengj commented Jan 20, 2017

For example, it would be good to have

B = A'
B * C, C * B, B * B

use BLAS if A and C are Matrix{Float64} or Matrix{Complex{Float64}}.

Oh, wait, I forgot that we don't have TransposedMatrix yet.

@andyferris
Copy link
Member Author

Right, the next step is a much more disruptive (in terms of code) PR that introduces TransposedMatrix.

Because of the combinatorial explosion of methods of having all the Ac_mul_B methods and dispatch on all the transpose and conjugation types, at that point I'll reverse tack and start removing any and all Ac_mul_B methods I see and re-write it in terms of a submethod of *(::ConjTransposedMatrix, ::AbstractMatrix) (or vectors, etc). (At that point it is natural that we look at deprecating the special parsing of a' * b.)

I'm guessing this is post v0.6 (primarily for the reason is that we've already run out of time) unless people really think it's worth pursuing now.

@andyferris
Copy link
Member Author

OK I've done a little bit more work on this and rebased... just the sparse module \ test to go, I think. That could be done separately or together, don't mind.

@inline similar(rowvec::RowVector) = RowVector(similar(parent(rowvec)))
@inline similar{T}(rowvec::RowVector, ::Type{T}) = RowVector(similar(parent(rowvec), transpose_type(T)))

# Resizing similar currently looses its RowVector property.
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

loses

@StefanKarpinski StefanKarpinski added this to the 0.6.0 milestone Jan 25, 2017
@StefanKarpinski
Copy link
Member

There are three reviewers, none of which have officially approved this. @stevengj, @tkelman, @Sacha0 – can you approve if this is good to go, modulo typos?

@andyferris
Copy link
Member Author

OK, I fixed the broken test (upgraded from "brittle" to "somewhat brittle") and typo, and have rebased.

Regarding the StridedArray interface functions (pointer and stride) - do people thin that these should be implemented here even though they won't be used, or should we fix/attack the StridedArray definition first (either by making it a trait or expanding its union)?

To me, this had felt like work that could be done along with TransposedMatrix - they will both be touching a similar set of function definitions (remember that RowVector simply unwraps itself for *, \, /, ..., but TransposedMatrix won't have that luxury if we get rid of Ac_mul_B, etc).

@inline ctranspose{T<:Real}(vec::AbstractVector{T}) = RowVector(vec)

@inline transpose(rowvec::RowVector) = rowvec.vec
@inline transpose(rowvec::ConjRowVector) = copy(rowvec.vec) # remove the ConjArray wrapper from any raw vector
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does this make a copy?

Copy link
Member Author

@andyferris andyferris Jan 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The logic was to make sure that vectors are always unwrapped, so that calls to BLAS are generally used. This is partly because the PR doesn't create a bunch of new methods for *(conjarray, array), etc. We also have "view semantic" for row vectors and "copy semantic" for vectors, which seemed more consistent than views for vectors and rowvectors and copies for everything else.

(That is, from the user's perspective, only the new RowVector type gets the view semantic, conjugated or not. The old types behave as before, so hopefully this is more backward compatible. Ideally these changes would have been made at the same time as TransposedMatrix to make all of this more consistent...)

Of course I open to doing it the other way, if people thing that makes more sense.

@@ -147,6 +154,11 @@ end

# Multiplication #

# inner product -> dot product specializations
@inline *{T<:Real}(rowvec::RowVector{T}, vec::AbstractVector{T}) = dot(parent(rowvec), vec)
@inline *(rowvec::ConjRowVector, vec::AbstractVector) = dot(rowvec', vec)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm worried about potential performance regressions if we don't also have *(rowvec,Matrix) = At_mul_B and *(conjrowvec,Matrix) = Ac_mul_B

Copy link
Member Author

@andyferris andyferris Jan 27, 2017

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The first is true, but not the way you think: you end up on this line so you get *(rowvec, matrix) = RowVector(At_mul_B(matrix, parent(rowvec))). Would that be a performance penalty?

I should chase down all the conjugated cases and make sure they do something sensible. I think what happens is a conjugated copy of the parent vector is created as a temporary, and the above gets called. (We should be able to remove this temporary - one way would have the output become a ConjRowVector which seems a little unusual).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

RowVector(At_mul_B(matrix, parent(rowvec))) should be fine.

@StefanKarpinski
Copy link
Member

What do we need to do to make forward progress here?

@andyferris
Copy link
Member Author

andyferris commented Jan 27, 2017

  • There might be a regression in vec' * mat such that an extra conjugated vector copy is made. Can I have feedback on whether people feel this will be a large or small problem, or could be addressed separately vs should be addressed in this PR? (This regression was introduced by Introduce RowVector as the transpose of a vector #19670)

  • Similarly, can we expect any performance difference between RowVector(At_mul_B(mat, vec)) and the old implementations of At_mul_B(vec, mat)? (thanks @stevengj)

  • Can I have a final word about providing the strided array functions (pointer and/or stride) for these views (given they wouldn't be used by code in Base at the moment)?

@stevengj
Copy link
Member

We can't have pointer etc. for ConjArray, and I think we should punt on that for RowVector (given that StridedArray is not yet a trait anyway).

end

@testset "RowVector conjugates" begin
v = [1+im, 1-im]
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

4-space

@andyferris
Copy link
Member Author

Rebased, indentation fixed.

(If you want to merge this now before "alpha" then we can look at any remaining Ac_mul_B performance regressions where a RowVector is returned which were introduced in #19670 in a separate, more focused PR. I don't see such performance penalties as being very large or very difficult to fix.)

@andyferris
Copy link
Member Author

Bump

@tkelman
Copy link
Contributor

tkelman commented Feb 2, 2017

I like it. If ConjMatrix isn't going to occur in normal operations yet, should we maybe hold off on exporting that alias for now?

edit: I also wish nanosoldier were working on master right now

@andyferris
Copy link
Member Author

I like it.

Thanks. I'm already toying with TransposedMatrix and refactoring matmul, etc. I think we'll see nice simplifications out of that.

If ConjMatrix isn't going to occur in normal operations yet, should we maybe hold off on exporting that alias for now?

Right. On one hand, I'm definitely already using it in my toying around, but it is one more symbol in the namespace. I kind-of felt having the vector and array forms and not the matrix form as being a bit odd and distinct to how these things typically come in triples. I'm easy.

@stevengj
Copy link
Member

stevengj commented Feb 2, 2017

LGTM.

I'd still like a method to improve performance in y = x'; x * A, but that can be a separate PR.

A lazy-view wrapper of an `AbstractArray`, taking the elementwise complex
conjugate. This type is usually constructed (and unwrapped) via the `conj()`
function (or related `ctranspose()`), but currently this is the default behavior
for `RowVector` only.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps a doctest here with an example of behaviour?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'll have to learn what a doctest is, first. I guess I'll try the docs.

"""
conj(rowvector)

Returns a `ConjArray` lazy view of the input, where each element is conjugated.
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Perhaps a @ref to ConjArray and maybe a doctest?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I didn't know about @ref either, but that seems useful.

@andyferris
Copy link
Member Author

andyferris commented Feb 5, 2017

Done, @KristofferC . I had trouble getting either of

[`'`](@ref)
[`.'`](@ref)

to work, though, but perhaps this is good enough.

@andyferris
Copy link
Member Author

I think this might be good to go?

@stevengj
Copy link
Member

stevengj commented Feb 6, 2017

Looks like the tests need to be updated after #20423.

@StefanKarpinski
Copy link
Member

Bump. This is one of the last remaining feature changes for 0.6.

@andyferris
Copy link
Member Author

@StefanKarpinski waiting on CI, but this should be good to merge now.

@@ -237,19 +235,3 @@ end

@test_throws DimensionMismatch ut\rv
end

# issue #20389
@testset "1 row/col vec*mat" begin
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

shouldn't delete these

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think one of the rebases must have done something odd here. I'll fix it.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

fixed

By default, this is used only in conjugation with `RowVector`, so that
both `transpose(vec)` and `ctranspose(vec)` both return views.
@andyferris
Copy link
Member Author

Squashed and rebased, for convenience.

@StefanKarpinski
Copy link
Member

So the pending work here is all adding optimized internal methods, right?

@andyferris
Copy link
Member Author

That's correct, Stefan. It shouldn't affect the user interface or behaviour.

@StefanKarpinski
Copy link
Member

Great – let's pull the trigger once CI passes everywhere.

@andyferris
Copy link
Member Author

Looks like Travis finally started.

BTW, Stefan, I have started on a branch here with TransposedMatrix and a refactoring (i.e. removal) of all the Ac_mul_B and Ac_mul_B! methods into * and A_mul_B!. So far I've attacked the generic and strided/BLAS multiplications, but that leaves symmetric, hermitian, diagonal, bidiagonal, triangular, sparse, etc, as well as the division methods... :) I've also made some movements of code so that BLAS is more modular. (I've noticed there are assumptions everywhere that the external libraries are present and it will take quite a bit of effort to turn BLAS, LAPACK, ARPACK, etc into proper modules that can be loaded at the end of LinAlg - ideally, optionally, with the caveat that there aren't native implementations of many things.)

# a "view" for conj(::AbstractArray)
@inline conj(rowvec::RowVector) = RowVector(conj(rowvec.vec))
"""
conj(rowvector)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

not worth going through CI again right now, but this should probably say conj(v::RowVector) to be specific?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I took the liberty of changing this over #20446

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@tkelman are we moving towards these more accurate type signatures in docs? I was copying the style that I've seen in a lot of other places, but for myself I would have more naturally written it like you have here.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

if the docstring description applies generally, then the signature should be general. if it's about a method for a specific type, then the signature should show that

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

thanks

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
linear algebra Linear algebra
Projects
None yet
Development

Successfully merging this pull request may close these issues.

8 participants