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

Restrict +,- and / to algebraic operations (for matrices). Use . versions for elementwise operations. Add UniformScaling matrix. Was:Don't let division with array default to elementwise division. #5810

Merged
merged 1 commit into from
Mar 10, 2014

Conversation

andreasnoack
Copy link
Member

Define Number rhs for Matrix lhs when dimensions match.

This is a followup on the discussion in #5807.

@jiahao
Copy link
Member

jiahao commented Feb 13, 2014

There is a test whose behavior needs to be changed too.

@stevengj
Copy link
Member

My inclination would be to not define Number/Matrix at all. Force the user to explicitly indicate either ./ or inv.

@goszlanyi
Copy link

Please don't do this.
And yes, force the user to be explicit about it.

@lindahua
Copy link
Contributor

Agree with not defining Number / Matrix. We should encourage user to make his meaning explicit when there may be two ways to understand the thing.

@andreasnoack
Copy link
Member Author

I'd be okay with removing Number/Matrix. The reason I included them was that we try to treat numbers a length one vectors in the rest of the linear algebra code. The example I used in #5807 was

julia> [1 1]\1
1x2 Array{Float64,2}:
 1.0  1.0

julia> [1 1]\[1]
2-element Array{Float64,1}:
 0.5
 0.5

@andreasnoack
Copy link
Member Author

I have tried to follow the suggestions in #5807. So now you can do

julia> A=randn(4,4);
julia> λ=eigvals(A)[1];
julia> det(A-λ*I)
5.0484003519825854e-15 + 0.0im

and

julia> zeros(3,3)+I
3x3 Array{Float64,2}:
 1.0  0.0  0.0
 0.0  1.0  0.0
 0.0  0.0  1.0

but get

julia> zeros(3,3)+1
WARNING: A::Array + x::Number is deprecated, use A .+ x instead.
 in + at deprecated.jl:19
3x3 Array{Float64,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

I have defined 1/A==inv(A) for now even though we were divided on the subject. As pointed out by at least @lindahua, this is a silent breakage in 1/A because it used to work elementwise. I still buy @toivoh's argument that (A*1)/A==A*(1/A), but admit the problem in silently changing the method. Maybe a solution could be to remove the method and then reintroduce it in its new definition, when people have had time to remove the old use of /.

During the adjustment of code in base and test I removed a few lines in mapslices which made its return type depend on the value of dims. Now mapslices(sum,Matrix,(1,2))==sum(Matrix,(1,2)). This change might be of interest to @JeffBezanson.

@lindahua
Copy link
Contributor

Any change of behavior should not be silent. The behavior of 1/A remains controversial, therefore a prudent approach should be to first deprecate this syntax for a period of time, and we may always reintroduce it using the algebraic semantics later if people agree to do so.

@@ -0,0 +1,49 @@
import Base: getindex, show, +, *
immutable IdentityMatrix{T<:Number} <: AbstractMatrix{T}
Copy link
Member

Choose a reason for hiding this comment

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

It still seems weird to me to call this IdentityMatrix because of the scale factor. Maybe ScaledIdentity or ScalingOperator?

Copy link
Member Author

Choose a reason for hiding this comment

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

Ah. I forgot to mention that. I remember you comment in #5807 and I agree that the name IdentityMatrix is not precise because of the λ, but couldn't come up with a better name. I think ScalingOperator/Matrix is used for a matrix that allow distinct diagonal elements. ScaledIdentity is fine.

Copy link
Contributor

Choose a reason for hiding this comment

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

+1 for ScaledIdentity.

Copy link
Contributor

Choose a reason for hiding this comment

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

@StefanKarpinski do you have any suggestion of this name?

Copy link
Member

Choose a reason for hiding this comment

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

ScaledIdentity seems reasonable.

Copy link
Contributor

Choose a reason for hiding this comment

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

What about Homothety ? Maybe it's just in french but it is quite common in math in my experience.

Copy link
Member

Choose a reason for hiding this comment

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

I've never heard that word before. Having looked it up, this isn't actually what it means – it's an affine transformation that this would be the non-translation part of. In other words, this is a homethety with zero translation. In any case, I think it best not to use that term.

Copy link
Contributor

Choose a reason for hiding this comment

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

You are right of course, but usually it is implied to be a linear homothety in a linear algebra context. Anyway, if it is not familiar in english it's not worth it at all - apart from being nice to the ear :-)

Copy link
Member Author

Choose a reason for hiding this comment

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

French wikipedia defines it without the translation, so we could call it Homothétie, but I am not sure it is the best PR strategy for the new type. English wikipedia suggests the name UniformScaling.

Copy link
Member

Choose a reason for hiding this comment

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

I like UniformScaling best so far. There's something a bit oxymoronic about ScaledIdentity, after all.

@andreasnoack
Copy link
Member Author

@stevengj Thank you for the comments. I have updated the pr.

Changes:

  • Deprecate Number/Matrix, added new version but commented out
  • UniformScaling instead of IdentityMatrix
  • Methods for \
  • Some documentation
  • Tests

import Base: getindex, showarray, +, -, *, /, ctranspose, transpose
import Base.LinAlg: SingularException
immutable UniformScaling{T<:Number} <: AbstractMatrix{T}
λ::T
Copy link
Contributor

Choose a reason for hiding this comment

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

Looks like this source file uses tab for indentation. I think the Julia base uses 4-spaces.

@toivoh
Copy link
Contributor

toivoh commented Mar 2, 2014

Btw, thanks for taking the lead on this!

@andreasnoack
Copy link
Member Author

Ok. I have changed matrix to operator in the documentation and tabs to spaces in uniformscaling.jl. Any further comments before merging?

@toivoh You're welcome. I agree that it would be very useful if the abstractions were documented. There has been several issues related to different perceptions of abstractions.

@lindahua
Copy link
Contributor

I think we should merge this as soon as possible, preferably before 0.3 release.

@JeffBezanson
Copy link
Member

+1

@JeffBezanson JeffBezanson added this to the 0.3 milestone Mar 10, 2014
@jiahao
Copy link
Member

jiahao commented Mar 10, 2014

Looks good to me.

…operator I. Let Number/Matrix be the inverse. Update docs.
andreasnoack added a commit that referenced this pull request Mar 10, 2014
Restrict +,- and / to algebraic operations (for matrices). Use . versions for elementwise operations. Add UniformScaling matrix.
@andreasnoack andreasnoack merged commit 4e1e0ab into master Mar 10, 2014
@andreasnoack andreasnoack deleted the anj/arraydiv branch March 10, 2014 06:23
@andreasnoack
Copy link
Member Author

I wasn't sure if this was a 0.3 thing, but now it is.

@@ -189,6 +189,16 @@ export PipeString
@deprecate svdfact(A,thin) svdfact(A,thin=thin)
@deprecate svdfact!(A,thin) svdfact(A,thin=thin)
@deprecate svd(A,thin) svd(A,thin=thin)
@deprecate (+)(A::Array{Bool},x::Bool) A .+ x
Copy link
Member

Choose a reason for hiding this comment

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

These deprecations will have to be replaced with a gentle error explaining why you can't do [1:10] + 4, when the deprecation period is over. A NoMethodError will be confusing and cause people to submit PRs with the "missing" methods.

Copy link
Member

Choose a reason for hiding this comment

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

@johnmyleswhite
Copy link
Member

I'm with @andreasnoackjensen on this one.

@StefanKarpinski
Copy link
Member

The place I've found it most annoying is when doing vector + scalar, where arguably, there is no ambiguity.

@nalimilan
Copy link
Member

Agreed. In that case you don't gain anything by forcing users to use .+ instead of +.

@BobPortmann
Copy link
Contributor

I don't like having to type .* to get elementwise multiplication. But as long as that is necessary I think the system should be consistent and .+ should be necessary for elementwise addition as well.

@stevengj
Copy link
Member

You don't have to type * to get elementwise multiplication for array*scalar.

(And to the purist argument that "linear algebra doesn't define matrix + scalar", I would say that (a) Julia has never been restricted to solely the operations of linear algebra, (b) math doesn't prevent you from defining whatever notation you want, and (c) the notation array+scalar has already been defined and is widespread in scientific computing. The only question is whether it is confusing, but the widespread adoption of array+scalar means that not defining it is far more confusing.)

@BobPortmann
Copy link
Contributor

Yes, but I thought before this change + worked elementwise between arrays but it didn't broadcast.

@lindahua
Copy link
Contributor

I don't have strong preferences either way.

If we restore the original behavior, then there would be some subtle differences between + and .+: the former supports broadcasting a scalar, while the latter supports broadcasting in a more flexible way (e.g, rand(3, 4) .+ rand(3)). Two broadcasting operators with such subtle different behaviors might introduce additional confusion, especially to beginners.

@stevengj
Copy link
Member

@BobPortmann, no, before this change, array+scalar was broadcasting, just like in Matlab.

@stevengj
Copy link
Member

@lindahua, couldn't you make exactly the same argument for array*scalar vs. array.*scalar? There will always be some overlap between the dot and non-dot operators.

The point is simply that +, like *, is not a general broadcasting operation: the only broadcasting-like operation it supports is array+scalar (along with addition of two arrays of the same shape). Considering that most beginners will be coming from languages where this is true, I find it hard to believe that requiring .+ will be less confusing.

@BobPortmann
Copy link
Contributor

@stevengj Yes, of course. I meant "array broadcasting" (i.e., modifying an array to match another array). "Scalar broadcasting" seems so natural I don't even think of it as broadcasting (although it is).

I'm fine with relaxing the rules for array + scalar, but really don't think it matters much. Once you get in the habit of typing all those annoying .'s you might as well keep typing them :-)

@goszlanyi
Copy link

@stevengj
I truly appreciate your momentum on this matter.
Two months ago (#6417) I thought that there was no hope.

@stevengj
Copy link
Member

@BobPortmann, as made concrete in #7226, the only controversy here is the behavior for array+scalar. No one is proposing that array+array work unless the arrays have the same shape.

@stevengj
Copy link
Member

@GaborOszlanyi, I think that a couple of months ago, many of us felt that .+ would become second nature after using it a few times (much like im vs. i or a[i] vs. a(i)).

Now that I've had a chance to use it, I'm sad to say that [1,2,3] .+ 4 has not grown on me. It's just a constant annoyance to be nagged by Julia to type a . when [1,2,3] + 4 is perfectly unambiguous. My sense from this discussion is that many other (but not all) Julia developers have come to the same conclusion.

Julia should be fun. Being nagged about pointless algebraic purity is not fun.

@StefanKarpinski
Copy link
Member

Language design, turns out to be a subfield of psychology, not mathematics – it's highly experimental – sometimes you've just got to try things out and see if they work or not. Even after a few months, this is still annoying.

@Keno
Copy link
Member

Keno commented Jun 11, 2014

Should we try enabling it for vectors only to see how that feels?

@StefanKarpinski
Copy link
Member

Seems too fiddly. I'm inclined to call this a failed experiment, re-allow scalar + array and be done with it.

@IainNZ
Copy link
Member

IainNZ commented Jun 12, 2014

With all due respect to theoretical purity of the distinction between 2, I prefer the old behavior, so I guess thats another 👍. I still forget it every time...

@ghost
Copy link

ghost commented Jun 12, 2014

Considering the amount of accidental pain that broadcasting in NumPy has caused me over the years I really like the current clear distinction in Julia. When working with arrays the only "type-safety" that you have are the dimensions of your arrays and whether or not something is a scalar, so I'll gladly take any additional safety that I can get.

@ViralBShah
Copy link
Member

I do like the clear separation between broadcasting and algebraic operations as @andreasnoackjensen and @johnmyleswhite suggest. I also like the fact that the accidental pain of broadcasting as in Numpy is avoided as @ninjin says. But, simultaneously, I really hate doing .+ for scalar + array.

I personally prefer that we only allow + for scalar + array, and not extend to other forms of broadcasting such as vector + array, which does feel wrong. At least, our scalars and arrays are cleanly separated, and hence we could get away with special behaviour for scalars. For all other forms of broadcasting, I do prefer .+.

@GunnarFarneback
Copy link
Contributor

Yes, please allow scalar + array. I use it all the time and for existing code the only feasible way was to override the deprecations and bring back scalar + array myself. For new code I've been using .+ and it has been a pain.

@JeffBezanson
Copy link
Member

It looks like the consensus leans in favor of bringing back scalar + array.

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

Successfully merging this pull request may close these issues.