re-allow vector*(1-row matrix) and transpose thereof #20423

stevengj · 2017-02-03T02:06:39Z

As discussed in JuliaLang/LinearAlgebra.jl#400, I think that forbidding this case will cause more confusion than it is worth.

StefanKarpinski · 2017-02-03T04:33:53Z

base/linalg/rowvector.jl

 *(::RowVector, ::RowVector) = throw(DimensionMismatch("Cannot multiply two transposed vectors"))
 @inline *(vec::AbstractVector, rowvec::RowVector) = vec .* rowvec
 *(vec::AbstractVector, rowvec::AbstractVector) = throw(DimensionMismatch("Cannot multiply two vectors"))
-*(mat::AbstractMatrix, rowvec::RowVector) = throw(DimensionMismatch("Cannot right-multiply matrix by transposed vector"))


Is this implemented somewhere?

So you want to allow this one too, Steven? In this case it would return a matrix, another generalized form of the outer product.

@StefanKarpinski, yes, it is implemented and I included a test: it falls back on the generic matrix×matrix since RowVector is a subtype of AbstractMatrix. I allow it because it is the transpose of the other case that we allow, and it is weird to allow one but not the other (x*A but not A'*x').

tkelman · 2017-02-03T05:41:20Z

We don't usually implicitly convert between different dimensionalities all that many places, I think we should deprecate this. Otherwise are we going to allow matrix * rowvector too, if the matrix happens to have 1 column?

andyferris · 2017-02-03T07:47:25Z

Tony, Stefan - to clarify, the line removed here will enable matrix * rowvector -> matrix.

(While I am in favour of deprecation, we do at least need to reimplement the method to go with that.)

stevengj · 2017-02-03T13:10:50Z

@tkelman, there was a long discussion about this in the other issue. Short answer (a) we have allowed this particular thing for a long time and many users expect it and (b) yes, this PR allows matrix×rowvector for 1-column matrices too — that is even more "natural" in terms of code — it requires no special methods, because rowvector is already a subtype of matrix. (Disabling matrix×rowvector requires special-case code.)

(Also (c): any short deprecation message will simply cause intense confusion for people coming from fields where they are used to thinking of vectors as being equivalent to 1-column matrices when doing linear algebra and (d) it's hard to properly say what to replace such code with without knowing where the 1-column matrix came from: e.g. x * ones(n,1) should be replaced with x * ones(n)'. A generic deprecation would probably call reshape, and that would give users the misleading impression that we are requiring them to write ugly reshape code for what should be a simple operation.)

tkelman · 2017-02-03T13:25:10Z

We haven't had a rowvector type for all this time that we've allowed dimensionally questionable operations. It's not that hard for people to grasp that in Julia a 1D vector and a 1 column matrix are different types and behave differently. I thought much of the point of 4774 is to allow us to make the distinction clearer.

stevengj · 2017-02-03T13:37:05Z

@tkelman, just because they are different types and can follow different rules doesn't mean that they should. Many people learned linear algebra by thinking of vectors as following the same rules as one-column matrices, and think of A*x in that way, and learn the outer product x*y' as following the matrix*matrix rules for 1col*1row. Why can't we support this, rather than forcing users to comprehend a higher level of abstraction in thinking of their vector spaces?

A RowVector is a subtype of AbstractMatrix that acts in other ways like 1-row matrix, and vector*rowvector is defined, so in that sense also it is weird for that to work but vector*(1-row matrix) to not work.

What do you gain by imposing your notion of purity on these operations? A lot of users will be confused and inconvenienced; what is the corresponding gain?

…matrices

tkelman · 2017-02-03T13:54:06Z

Using the type system to write code that more clearly reflects when it does and does not work, and catches bugs earlier by allowing fewer questionable operations. Users having to think about the dimensionality of their arrays isn't a big ask. The vector-vs-matrix distinction in Julia is one of the many useful differences between us and Matlab, and not all that abstract or complicated.

stevengj · 2017-02-03T13:55:26Z

(Another clue: this PR is a net deletion of code, even if you count the additional tests. Supporting fewer operations here requires more code, which is another clue that a restrictive view of matvec here isn't worth it.)

stevengj · 2017-02-03T13:56:41Z

@tkelman, this is not a Matlab thing. This is a way that some people think about linear algebra. Why do we need to try (and probably fail) to convince them that their way of thinking about linear algebra is "wrong"? (They will no doubt conclude instead that Julia is overly pedantic and annoying.)

tkelman · 2017-02-03T13:59:24Z

In terms of realizing that linear algebra as code, matlab has no distinction, we do. The question is whether we want that distinction to be reflected in what operations we allow, or partially try to hide it only for matvec.

stevengj · 2017-02-03T14:03:30Z

Complex and real numbers are different types in Julia. Are we "partially hiding" that distinction by allowing complex*real? I don't think this is the right way to frame the question.

Many people think of linear algebra as allowing them to multiply vectors and one-row matrices. It doesn't cost us anything to support that viewpoint, and in fact it saves us code. It will confuse them if we forbid it, not because of Matlab, but because this is how they think of matvec operations.

stevengj · 2017-02-03T14:55:11Z

See also JuliaLang/LinearAlgebra.jl#400 by @StefanKarpinski

fredrikekre · 2017-02-03T17:14:47Z

Should it be more than twice as slow? That is also something to consider for people that would use this method.

julia> using BenchmarkTools

julia> for N in (10, 100, 1_000, 10_000)
           times = zeros(2, length(N))
           x = rand(N); y = rand(N);

           yt = y'; ymat = reshape(y, 1, length(y));

           x*yt;
           x*ymat;

           @btime $x*$yt;
           @btime $x*$ymat;
       end
  163.458 ns (1 allocation: 896 bytes)
  326.870 ns (3 allocations: 992 bytes)
  5.138 μs (2 allocations: 78.20 KiB)
  9.544 μs (4 allocations: 78.30 KiB)
  556.765 μs (2 allocations: 7.63 MiB)
  1.333 ms (4 allocations: 7.63 MiB)
  88.486 ms (2 allocations: 762.94 MiB)
  289.008 ms (4 allocations: 762.94 MiB)

mbauman · 2017-02-03T17:39:41Z

I'm in favor of this PR for now. We should re-evaluate this in 0.6 to make a consistent decision alongside trailing singleton dimensions and partial indexing. Either our we allow our arrays to masquerade as other dimensionalities as convenient or we don't.

stevengj · 2017-02-03T20:00:41Z

@fredrikekre, the speed is a separate issue; it's not intrinsic to what operations we allow (since the underlying calculations are essentially identical), it's just a question of what specialized methods we provide.

(Right now, if you treat a RowVector as a matrix in a matmul, it goes through generic code that is not as optimized. I think it is fine for the less-preferred way of expressing something to be slower by a constant factor ... the kinds of operations we're talking about are often not that performance critical anyway.)

StefanKarpinski · 2017-02-03T22:56:59Z

We should re-evaluate this in 0.6 to make a consistent decision alongside trailing singleton dimensions and partial indexing. Either our we allow our arrays to masquerade as other dimensionalities as convenient or we don't.

Note that allowing indexing with 1 into non-existent dimensions is quite different than partial linear indexing – in particular, partial linear indexing leaks (possibly invalid) assumptions about memory layout of arrays, whereas allowing indexing into non-existing dimensions with a 1 does not.

mbauman · 2017-02-03T23:17:20Z

Yes, but trailing singleton dimensions are quite alike inserting implicit singleton dimensions when you use more than one but fewer than N indices to index into an N-dimensional array. That's what I meant by partial indexing. We need better vocabulary here.

StefanKarpinski · 2017-02-03T23:21:12Z

Ok, let's go with this now and then reconsider the whole thing in the 0.6 cycle.

andyferris · 2017-02-04T00:04:31Z

OK, I agree to merge this for v0.6 and revisit trailing dimensions and this together in the next development cycle (I think implied trailing singleton dimensions may be a hole in the argument in this comment).

@stevengj you missed some methods in other files (triangular and diagonal - see #20429 for exactly where), which really should be dealt with in this PR. There may also be other matrix * rowvector disambiguity methods splattered around linalg/ and sparse/, not sure right now (sorry I don't have much time at the moment).

Complex and real numbers are different types in Julia. Are we "partially hiding" that distinction by allowing complex*real? I don't think this is the right way to frame the question.

Here, the interface is defined by Number. The question is always which types (or traits) define the interface, and what methods that interface should consist of, and how does this implement concepts users need? We are talking about what interfaces are supported by AbstractVector, AbstractMatrix and RowVector.

…lar cases

stevengj · 2017-02-04T01:45:10Z

@andyferris, done. Note that since diagonal and triangular matrices in Julia are always square, the only function of those methods was to throw an error in the 1x1 case, which hardly seems worth the trouble.

tkelman · 2017-02-04T09:39:14Z

base/linalg/matmul.jl

@@ -80,6 +80,7 @@ function (*){T,S}(A::AbstractMatrix{T}, x::AbstractVector{S})
    TS = promote_op(matprod, T, S)
    A_mul_B!(similar(x,TS,size(A,1)),A,x)
 end
+(*)(A::AbstractVector, B::AbstractMatrix) = reshape(A,length(A),1)*B


so this method didn't previously exist?

It did. AFAICT, it had a similar implementation modulo inlining. (I have trouble with the mobile version of GitHub, sorry).

andyferris · 2017-02-04T10:18:02Z

It did. AFAICT it had a similar implementation modulo inlining.

andyferris · 2017-02-04T10:22:51Z

@stevengj Thanks, it's still good to remove those error messages I wrote, since they would be non-sensical.

andyferris · 2017-02-04T13:29:07Z

@stevengj I think that the tests can be fixed by similar test changes as in #20429.

stevengj · 2017-02-04T16:08:54Z

@andyferris, I think it should be fixed now—I fixed it by adding generic methods, rather than specialized methods for triangular, and deleting some specialized triangular methods.

andyferris · 2017-02-04T21:22:55Z

base/linalg/matmul.jl

@@ -81,6 +81,13 @@ function (*){T,S}(A::AbstractMatrix{T}, x::AbstractVector{S})
    A_mul_B!(similar(x,TS,size(A,1)),A,x)
 end

+# these will throw a DimensionMismatch unless B has 1 row (or 1 col for transposed case):
+A_mul_B(a::AbstractVector, B::AbstractMatrix) = A_mul_B(reshape(a,length(a),1),B)


A_mul_B? Is this meant to be *?

whoops, right, this doesn't exist

andyferris · 2017-02-04T21:25:15Z

base/linalg/matmul.jl

+A_mul_B(a::AbstractVector, B::AbstractMatrix) = A_mul_B(reshape(a,length(a),1),B)
+A_mul_Bt(a::AbstractVector, B::AbstractMatrix) = A_mul_Bt(reshape(a,length(a),1),B)
+A_mul_Bt(A::AbstractMatrix, b::AbstractVector) = A_mul_Bt(A,reshape(b,length(b),1))
+A_mul_Bc(a::AbstractVector, B::AbstractMatrix) = A_mul_Bc(reshape(a,length(a),1),B)


Should At_mul_B(::RowVector, ::AbstractMatrix) and friends be part of this list too? Or I guess it works as generic matrix-matrix multiplication?

Since RowVector is an AbstractMatrix, it should always work with the generic fallback.

andyferris · 2017-02-04T22:30:02Z

Cool - now LGTM. Thanks Steven.

re-allow vector*(1-row matrix) and transpose (closes #20389)

7f977ed

stevengj added the linear algebra Linear algebra label Feb 3, 2017

StefanKarpinski reviewed Feb 3, 2017

View reviewed changes

This was referenced Feb 3, 2017

Include deprecated methods for vector * matrix #20429

Closed

0.6: vec*mat throws "Cannot left-multiply a matrix by a vector" even when mat is 1 x n JuliaLang/LinearAlgebra.jl#400

Closed

delete more methods that are just there to throw errors on 1 row/col …

30f9cd1

…matrices

rm more DimensionMismatch methods for RowVector: diagonal and triangu…

9e24fbf

…lar cases

tkelman reviewed Feb 4, 2017

View reviewed changes

add missing vector*matrix fallbacks

b799c16

andyferris reviewed Feb 4, 2017

View reviewed changes

A_mul_B does not exist, added missing A_mul_Bc case

c7ad5e6

stevengj merged commit e3bb870 into JuliaLang:master Feb 5, 2017

stevengj deleted the vecmat branch February 5, 2017 23:00

stevengj mentioned this pull request Feb 6, 2017

Added ConjArray wrapper type for conjugate array views #20047

Merged

andyferris pushed a commit to andyferris/julia that referenced this pull request Feb 8, 2017

Fix tests from change in JuliaLang#20423

fc962c5

andyferris mentioned this pull request Feb 16, 2017

A few mini tests for linalg #20613

Merged

mbauman mentioned this pull request Mar 28, 2019

Add some (matrix) functions for UniformScaling #28872

Merged

13 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

re-allow vector*(1-row matrix) and transpose thereof #20423

re-allow vector*(1-row matrix) and transpose thereof #20423

stevengj commented Feb 3, 2017 •

edited

Loading

StefanKarpinski Feb 3, 2017

andyferris Feb 3, 2017

stevengj Feb 3, 2017 •

edited

Loading

tkelman commented Feb 3, 2017

andyferris commented Feb 3, 2017

stevengj commented Feb 3, 2017 •

edited

Loading

tkelman commented Feb 3, 2017

stevengj commented Feb 3, 2017 •

edited

Loading

tkelman commented Feb 3, 2017

stevengj commented Feb 3, 2017

stevengj commented Feb 3, 2017 •

edited

Loading

tkelman commented Feb 3, 2017

stevengj commented Feb 3, 2017 •

edited

Loading

stevengj commented Feb 3, 2017

fredrikekre commented Feb 3, 2017

mbauman commented Feb 3, 2017 •

edited

Loading

stevengj commented Feb 3, 2017 •

edited

Loading

StefanKarpinski commented Feb 3, 2017

mbauman commented Feb 3, 2017

StefanKarpinski commented Feb 3, 2017

andyferris commented Feb 4, 2017 •

edited

Loading

stevengj commented Feb 4, 2017

tkelman Feb 4, 2017

andyferris Feb 4, 2017

andyferris commented Feb 4, 2017

andyferris commented Feb 4, 2017

andyferris commented Feb 4, 2017

stevengj commented Feb 4, 2017 •

edited

Loading

andyferris Feb 4, 2017 •

edited

Loading

stevengj Feb 4, 2017

andyferris Feb 4, 2017

stevengj Feb 4, 2017

andyferris commented Feb 4, 2017

re-allow vector*(1-row matrix) and transpose thereof #20423

re-allow vector*(1-row matrix) and transpose thereof #20423

Conversation

stevengj commented Feb 3, 2017 • edited Loading

StefanKarpinski Feb 3, 2017

Choose a reason for hiding this comment

andyferris Feb 3, 2017

Choose a reason for hiding this comment

stevengj Feb 3, 2017 • edited Loading

Choose a reason for hiding this comment

tkelman commented Feb 3, 2017

andyferris commented Feb 3, 2017

stevengj commented Feb 3, 2017 • edited Loading

tkelman commented Feb 3, 2017

stevengj commented Feb 3, 2017 • edited Loading

tkelman commented Feb 3, 2017

stevengj commented Feb 3, 2017

stevengj commented Feb 3, 2017 • edited Loading

tkelman commented Feb 3, 2017

stevengj commented Feb 3, 2017 • edited Loading

stevengj commented Feb 3, 2017

fredrikekre commented Feb 3, 2017

mbauman commented Feb 3, 2017 • edited Loading

stevengj commented Feb 3, 2017 • edited Loading

StefanKarpinski commented Feb 3, 2017

mbauman commented Feb 3, 2017

StefanKarpinski commented Feb 3, 2017

andyferris commented Feb 4, 2017 • edited Loading

stevengj commented Feb 4, 2017

tkelman Feb 4, 2017

Choose a reason for hiding this comment

andyferris Feb 4, 2017

Choose a reason for hiding this comment

andyferris commented Feb 4, 2017

andyferris commented Feb 4, 2017

andyferris commented Feb 4, 2017

stevengj commented Feb 4, 2017 • edited Loading

andyferris Feb 4, 2017 • edited Loading

Choose a reason for hiding this comment

stevengj Feb 4, 2017

Choose a reason for hiding this comment

andyferris Feb 4, 2017

Choose a reason for hiding this comment

stevengj Feb 4, 2017

Choose a reason for hiding this comment

andyferris commented Feb 4, 2017

stevengj commented Feb 3, 2017 •

edited

Loading

stevengj Feb 3, 2017 •

edited

Loading

stevengj commented Feb 3, 2017 •

edited

Loading

stevengj commented Feb 3, 2017 •

edited

Loading

stevengj commented Feb 3, 2017 •

edited

Loading

stevengj commented Feb 3, 2017 •

edited

Loading

mbauman commented Feb 3, 2017 •

edited

Loading

stevengj commented Feb 3, 2017 •

edited

Loading

andyferris commented Feb 4, 2017 •

edited

Loading

stevengj commented Feb 4, 2017 •

edited

Loading

andyferris Feb 4, 2017 •

edited

Loading