Make behavior of `^(::Matrix{<:Integer},::Integer)` consistent with behavior of `^(::Integer, ::Integer)` #23368

afniedermayer · 2017-08-20T15:26:23Z

Fixes JuliaLang/LinearAlgebra.jl#459. Makes ^(::Matrix{S},::T) where {S<:Integer, T<:Integer} type stable. Also ensures the same promotion rules as ^(::S, ::T).

… Also enure that the same promotion rules as `^(::S, ::T)`

…ayer/julia into afniedermayer/fix23366

…atrices

andreasnoack · 2017-08-21T02:12:09Z

Thanks for opening the PR. I think it would be a good idea simply to make power_by_squaring work for matrices. Then the matrix function could just call power_by_squaring directly when the power is an integer. What needs to be fixed is just the tests in

julia/base/intfuncs.jl

Lines 177 to 178 in 7a1156c

    
           x == 1 && return copy(x) 
        
           x == -1 && return iseven(p) ? one(x) : copy(x)

. You can use isone instead of looping as you do here. I think it is fine to use isone also when testing the negative identity. It is a corner case.

afniedermayer · 2017-08-21T06:01:35Z

@andreasnoack , thanks for reviewing this PR. Sure, simply using isone(A) || isone(-A) could be done, but it would come with a major performance penalty since you have to make a copy of the matrix -A. Consider this example (where the function _isabsone essentially does the same as the for loop I've written):

function _isabsone(A) # checks whether the absolute value of A isone
    m, n = size(A)
    if m != n || !isdiag(A)
        return false
    end
    is_eye, is_minus_eye = true, true
    for i = 1:m
        is_eye = is_eye && A[i, i] == 1
        is_minus_eye = is_minus_eye && A[i, i] == -1
        if !is_eye && !is_minus_eye
            return false
        end
    end
    return true
end
_isone_twice(A) = isone(A) || isone(-A)

A = fill(42,(10000,10000));
_isabsone(A); _isone_twice(A); # warm up JIT
@time _isabsone(A)
@time _isone_twice(A)

which gives the output:

  0.000004 seconds (4 allocations: 160 bytes)
  0.316946 seconds (6 allocations: 762.940 MiB, 1.18% gc time)

So _isone_twice is about 10000 times slower than _isoneabs and uses about 1GB more memory in this particular example.

If you think the simplicity is worth this performance penalty, I'll go ahead and change the pull request as you suggested.

fredrikekre · 2017-08-21T06:45:59Z

On another note, is it actually useful that we allow ^(x::Integer, p::Integer) when p < 0 if x == 1 or x == -1? If we throw here for all negative p we don't need to worry about checking isone(A) or isone(-A) for and integer matrix A :)

codecov-io · 2017-08-21T08:41:41Z

Codecov Report

Merging #23368 into master will decrease coverage by 0.02%.
The diff coverage is 66.66%.

@@            Coverage Diff             @@
##           master   JuliaLang/julia#23368      +/-   ##
==========================================
- Coverage   55.01%   54.98%   -0.03%     
==========================================
  Files         275      275              
  Lines       51145    51146       +1     
==========================================
- Hits        28138    28124      -14     
- Misses      23007    23022      +15

Impacted Files	Coverage Δ
base/path.jl	`55.46% <ø> (ø)`	⬆️
base/strings/util.jl	`64.33% <ø> (ø)`	⬆️
base/strings/utf8proc.jl	`55.17% <ø> (ø)`	⬆️
base/mpfr.jl	`51.98% <66.66%> (+0.11%)`	⬆️
base/distributed/messages.jl	`24.63% <0%> (-17.4%)`	⬇️
base/gcutils.jl	`14.28% <0%> (-7.15%)`	⬇️
base/distributed/remotecall.jl	`49.23% <0%> (-1.54%)`	⬇️
base/printf.jl	`60.3% <0%> (+0.15%)`	⬆️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 953bfc3...8267b1b. Read the comment docs.

afniedermayer · 2017-08-21T09:53:51Z

One of the tests fails on i686, but this seems to be the same test as reported in #23371. (I spent some time trying to figure out why the test fails before finding #23371. I wish @inferred gave more informative error messages...)

andreasnoack · 2017-08-22T03:00:35Z

I think it is worth to avoid the code duplication. With this change, the existing power_by_squaring is sufficient. The slowdown will only de significant for negative integer powers of the the negative identity. It is highly unlikely that anybody will need that to be fast.

…er},::Integer)` as for `^(::Integer,::Integer)`

afniedermayer · 2017-08-22T09:17:15Z

@andreasnoack , yes, I guess p<0 doesn't happen that often. I adjusted the code to avoid code duplication. If p<0 should turn out to be more common, one could still modify the code at some point in the future. (Or people could just use UniformScaling instead of a Matrix.)

fredrikekre · 2017-08-22T09:24:38Z

test/linalg/dense.jl

+                                            UInt128, UInt16, UInt32, UInt64, UInt8,
+                                            BigInt)
+        @test_throws DomainError elty[1 1;1 0]^-2
+        @test elty[1 1;1 0]^2 == elty[2 1;1 1]


Perhaps use some @inferred here since thats what this PR is about? Something like

@test (@inferred elty[1 1; 1 0]^2) == elty[2 1; 1 1]

should do the trick :)

Yes, you're right, I'll adjust that. (I'm not a big fan of @inferred in unit tests, since it results in confusing test failure messages, see https://discourse.julialang.org/t/more-informative-inferred-especially-for-unit-testing/5481 . But I guess there's no better solution at the moment.)

fredrikekre · 2017-08-22T09:26:58Z

test/linalg/dense.jl

+        end
+        for elty2 = (Int64, BigInt)
+            TT = Base.promote_op(^, elty, elty2)
+            @test Base.return_types(^, (Array{elty,2}, elty2)) == [Array{TT,2}]


I think this could also be rewritten using @inferred. (But this test might be redundant if you add an @inferred above)

The test isn't redundant because there's an additional issue going on that I discovered while fixing this issue. So far [1 1;1 0]^big(10000) returned Matrix{Int64} (and negative values) rather than promoting to Matrix{BigInt}. This was inconsistent with 2^big(10000). I fixed this issue as well.

Okay nice, perhaps worth adding a comment about that so it does not get removed by someone later on. Perhaps this could otherwise be tested with something like

@test ([1 1; 1 0]^big(1))::Matrix{BigInt} == [1 1; 1 0]

such that we test the return value and not just the type? :)

I think the test would only check for one particular value, whereas return_type proves type correctness for all possible values of a type. And [1 1;1 0]^big(1) is just one example of type promotion, there's others as well, e.g. Int32[1 1;1 0]^Int64(1). So it should be:

@test (@inferred elty[1 1;1 0]^elty2(1))::Matrix{TT} == [1 1;1 0]

I'll also add a comment to clarify.

fredrikekre · 2017-08-22T09:29:46Z

base/intfuncs.jl

-        x == 1 && return copy(x)
-        x == -1 && return iseven(p) ? one(x) : copy(x)
+        isone(x) && return copy(x)
+        isone(-x) && return iseven(p) ? one(x) : copy(x)
        throw_domerr_powbysq(p)


The error message printed in throw_domerr_powbysq could perhaps be updated to apply to integer matrices as well? Easiest would probably be to add a second method that explains the situation for integer matrices. Then this line could instead be

throw_domerr_powbysq(x, p)

and dispatch on x.

I had the same thought. I was a bit reluctant to make the error handling more complicated, since this A^p for p<0 seems to be a are thing for matrices. Where should throw_domerr_powbysq(x::AbstractMatrix,p) go? I guess dense.jl would be the right place, since intfucs.jl doesn't know anything about matrices.

I would put it here, since this is where it is used. If we add a specialized power by squaring method for matrices it can be moved to dense.jl later on. (Which, btw would be a good idea, since we can probably allocate some work arrays and allocate less in matrix multiplication, but thats for another PR :)).

afniedermayer · 2017-08-22T09:31:38Z

@fredrikekre , that's an important question you're asking: maybe p<0 should lead to a domain error even if isone(x) or isone(-x) both for AbstractArray{<:Integer}s and Integers. However, that would be a more radical change that might break existing integer code. So I think such a change would be more suitable for a separate PR. This PR just makes sure that AbstractArray{<:Integer} and Integer behave the same way.

`^(::Matrix{<:Integer},::Integer)`, added `@inferred` and comments to test

afniedermayer · 2017-08-22T11:55:26Z

There was a merge conflict because throw_domerr_powbysq(::Int64) had to be updated to throw_domerr_powbysq(::Int64, ::Int64) in the faq and the line numbers changed at the same time. I reverted my changes in the faq. The doctests for the faq would have to be fixed afterwards.

`^(::Matrix{<:Integer},::Integer)`, added `@inferred` and comments to test

…ayer/julia into afniedermayer/fix23366

andreasnoack · 2017-08-22T14:00:38Z

base/intfuncs.jl

-to_power_type(x) = x
-@noinline throw_domerr_powbysq(p) = throw(DomainError(p,
+to_power_type(x) = convert(Base.promote_op(*, typeof(x), typeof(x)), x)
+@noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p,
    string("Cannot raise an integer x to a negative power ", p, '.',


Since this is generic function, it might make sense to make this message generic if it doesn't make it (much) less informative. What about just writing "Cannot raise an integer x to a negative power. Try converting input to float". We could keep this version but restrict the inputs to Integers.

So you mean three versions of throw_domerr_powbysq: for ::Any (short message), ::Integer (message above), and ::AbstractMatrix (message below)?

andreasnoack · 2017-08-22T14:02:53Z

base/linalg/dense.jl

@@ -358,6 +358,10 @@ kron(a::AbstractVector, b::AbstractMatrix) = kron(reshape(a, length(a), 1), b)

 # Matrix power
 (^)(A::AbstractMatrix, p::Integer) = p < 0 ? Base.power_by_squaring(inv(A), -p) : Base.power_by_squaring(A, p)
+function (^)(A::AbstractMatrix{T}, p::Integer) where T<:Integer
+    TT = Base.promote_op(^, T, typeof(p))


Why is it necessary to promote here? I don't think this version is required. Couldn't you just change the previous method to

(^)(A::AbstractMatrix, p::Integer) = Base.power_by_squaring(A, p)

?

I guess the previous discussion got hidden by github: #23368 (comment)
The conversion is needed to make sure that Integer promotion rules carry over to matrices. For example, without this conversion [1 1;1 0]^big(10000) would return Matrix{Int64} (and negative values) rather than Matrix{BigInt}. This would be inconsistent with the behavior of 2^big(10000), which returns BigInt.

I see.

Below you can just write Base.power_by_squaring(convert(AbstractMatrix{TT}, A), p). Then it should be ready to go.

…sion

fredrikekre

Since this is breaking it would be nice to add a sentence in NEWS.md about this. It would be great if you could do that @afniedermayer, otherwise I (or someone else) can fix it in a separate PR.

afniedermayer · 2017-08-23T04:41:13Z

I added some text to NEWS.md. I also changed the title of this PR and of JuliaLang/LinearAlgebra.jl#459 to reflect that this is not just about type stability but also about promotion rules (e.g. [1 1;1 0]^big(2)).

We're good to go from my side. @fredrikekre and @andreasnoack , thanks for reviewing this PR!

fredrikekre · 2017-08-23T09:25:29Z

I split the NEWS.md bullet into two, and removed some unnecessary Base. qualifications, should be good to go when CI finishes 👍

fredrikekre · 2017-08-23T23:00:46Z

@andreasnoack wanna take a last look/merge?

afniedermayer added 6 commits August 20, 2017 16:58

Make ^(::Matrix{S},::T) where {S<:Integer, T<:Integer} type stable.…

119bb9d

… Also enure that the same promotion rules as `^(::S, ::T)`

Make ^(::Matrix{S},::T) where {S<:Integer, T<:Integer} type stable.…

b21baf2

… Also enure that the same promotion rules as `^(::S, ::T)`

Merge branch 'afniedermayer/fix23366' of https://github.com/afniederm…

c7d3de1

…ayer/julia into afniedermayer/fix23366

fix typo

c4e79d7

fixed typo

9fbe23f

Adjust tests for ^(::Matrix{Int64},Int64) for symmetric/hermitian m…

8267b1b

…atrices

kshyatt added linear algebra Linear algebra types and dispatch Types, subtyping and method dispatch labels Aug 20, 2017

avoid code duplication by using the same code for `^(::Matrix{<:Integ…

b5b3bc7

…er},::Integer)` as for `^(::Integer,::Integer)`

fredrikekre reviewed Aug 22, 2017

View reviewed changes

fredrikekre added breaking This change will break code needs news A NEWS entry is required for this change labels Aug 22, 2017

added specialized domain error message for

cd7a0a3

`^(::Matrix{<:Integer},::Integer)`, added `@inferred` and comments to test

afniedermayer added 3 commits August 22, 2017 13:58

added specialized domain error message for

56f5091

`^(::Matrix{<:Integer},::Integer)`, added `@inferred` and comments to test

Merge branch 'afniedermayer/fix23366' of https://github.com/afniederm…

ba5b3ec

…ayer/julia into afniedermayer/fix23366

reverted change in faq

76a33e2

andreasnoack reviewed Aug 22, 2017

View reviewed changes

different throw_domerr_powbysq for Any and Integer, simplified conver…

1a22961

…sion

fredrikekre approved these changes Aug 22, 2017

View reviewed changes

afniedermayer mentioned this pull request Aug 23, 2017

Behavior of ^(::Matrix{<:Integer}, ::Integer) inconsistent with behavior of ^(::Integer, ::Integer) JuliaLang/LinearAlgebra.jl#459

Closed

Adjusted "Breaking changes" section of news

e996b78

afniedermayer changed the title ~~Make ^(::Matrix{<:Integer},::Integer) type stable~~ Make behavior of ^(::Matrix{<:Integer},::Integer) consistent with behavior of ^(::Integer, ::Integer) Aug 23, 2017

fredrikekre added 3 commits August 23, 2017 11:18

split news item into two bullets

0bf799c

remove unnecessary Base. qualification

7fb880b

remove more unnecessary Base. qualifications

e1c334d

fredrikekre removed the needs news A NEWS entry is required for this change label Aug 23, 2017

andreasnoack merged commit 6378d7d into JuliaLang:master Aug 24, 2017

afniedermayer deleted the afniedermayer/fix23366 branch August 24, 2017 05:39

JeffBezanson removed the types and dispatch Types, subtyping and method dispatch label Oct 3, 2017

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Make behavior of `^(::Matrix{<:Integer},::Integer)` consistent with behavior of `^(::Integer, ::Integer)` #23368

Make behavior of `^(::Matrix{<:Integer},::Integer)` consistent with behavior of `^(::Integer, ::Integer)` #23368

afniedermayer commented Aug 20, 2017

andreasnoack commented Aug 21, 2017 •

edited by tkelman

Loading

afniedermayer commented Aug 21, 2017

fredrikekre commented Aug 21, 2017

codecov-io commented Aug 21, 2017

afniedermayer commented Aug 21, 2017

andreasnoack commented Aug 22, 2017

afniedermayer commented Aug 22, 2017

fredrikekre Aug 22, 2017

afniedermayer Aug 22, 2017

fredrikekre Aug 22, 2017

afniedermayer Aug 22, 2017

fredrikekre Aug 22, 2017

afniedermayer Aug 22, 2017

fredrikekre Aug 22, 2017

afniedermayer Aug 22, 2017

fredrikekre Aug 22, 2017

afniedermayer commented Aug 22, 2017

afniedermayer commented Aug 22, 2017 •

edited

Loading

andreasnoack Aug 22, 2017

afniedermayer Aug 22, 2017

andreasnoack Aug 22, 2017

afniedermayer Aug 22, 2017

andreasnoack Aug 22, 2017

fredrikekre left a comment

afniedermayer commented Aug 23, 2017

fredrikekre commented Aug 23, 2017

fredrikekre commented Aug 23, 2017

Make behavior of ^(::Matrix{<:Integer},::Integer) consistent with behavior of ^(::Integer, ::Integer) #23368

Make behavior of ^(::Matrix{<:Integer},::Integer) consistent with behavior of ^(::Integer, ::Integer) #23368

Conversation

afniedermayer commented Aug 20, 2017

andreasnoack commented Aug 21, 2017 • edited by tkelman Loading

afniedermayer commented Aug 21, 2017

fredrikekre commented Aug 21, 2017

codecov-io commented Aug 21, 2017

Codecov Report

afniedermayer commented Aug 21, 2017

andreasnoack commented Aug 22, 2017

afniedermayer commented Aug 22, 2017

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

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

Choose a reason for hiding this comment

afniedermayer commented Aug 22, 2017

afniedermayer commented Aug 22, 2017 • edited Loading

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

Choose a reason for hiding this comment

fredrikekre left a comment

Choose a reason for hiding this comment

afniedermayer commented Aug 23, 2017

fredrikekre commented Aug 23, 2017

fredrikekre commented Aug 23, 2017

Make behavior of `^(::Matrix{<:Integer},::Integer)` consistent with behavior of `^(::Integer, ::Integer)` #23368

Make behavior of `^(::Matrix{<:Integer},::Integer)` consistent with behavior of `^(::Integer, ::Integer)` #23368

andreasnoack commented Aug 21, 2017 •

edited by tkelman

Loading

afniedermayer commented Aug 22, 2017 •

edited

Loading