Add scale, scale! and copy! for triangular matrices #12683

mfasi · 2015-08-18T17:41:30Z

Add the methods scale! and copy! for triangular matrices. This addition is motivated by a few improvements suggested by @stevengj in #12584.

I tried to add the scale and copy methods as well. Not sure I did something reasonable there, as somewhere along the line I started listening to the compilation warnings without really understanding what was going on.

stevengj · 2015-08-18T17:51:22Z

base/linalg/triangular.jl

@@ -283,6 +283,155 @@ function -(A::UnitUpperTriangular)
    UpperTriangular(Anew)
 end

+# copy and scale
+function copy!(A::UpperTriangular,B::UpperTriangular)
+    copy!(A.data,B.data)


According to my tests, it is faster to just loop and copy the upper triangle.

mfasi · 2015-08-21T11:57:08Z

@stevengj Thank you for reviewing the patch. I tried my best to overcome the issues you pointed out. Two of them I don't know how to solve are:

The Real scaled by Complex version of the methods: without them, I get compilation warnings saying that my methods are ambiguous with

function scale{R<:Real,S<:Complex}(X::AbstractArray{R}, s::S)
    Y = Array(promote_type(R,S), size(X))
    copy!(Y, X)
    scale!(Y, s)
end

in base/linalg/generic.jl. I guess there must be a reason why this method is there and might not be wise to change it. Maybe the ambiguity can be solved in a somewhat cleverer manner, I don't really know.

The issue with BigFloat, on the other hand, depends -I think- on a feature of the Array(BigFloat) function. Indeed I have that

julia> Array(BigFloat, 5, 5)
5x5 Array{Base.MPFR.BigFloat,2}:
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef

and thus that

julia> A = convert(Matrix{BigFloat},randn(5,5)); similar(A)
5x5 Array{Base.MPFR.BigFloat,2}:
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef
 #undef  #undef  #undef  #undef  #undef

As the scale functions fill just the triangle, the remaining half-or-so of the matrix still has #undefs and the checks copy(A) == A fails, for example. I don't know how to handle this issue. A solution might be redefining the equality operator == for triangular matrices, but I'm not sure that would be completely safe.

Any suggestions?

andreasnoack · 2015-08-21T11:59:05Z

base/linalg/triangular.jl

+        end
+    end
+end
+function copy{T<:AbstractTriangular}(A::T)


Why is this method necessary? I think the copy method for AbstractArray can handle this one.

I hadn't noticed that, thank you!

mfasi · 2015-08-21T12:26:47Z

@andreasnoack Thank you for suggesting to fall back to scale and copy. I'm not sure this can handle the BigFloatissue though (see the inline comments).

andreasnoack · 2015-08-21T12:34:23Z

Sorry. The other comments had been folded because of new pushes to the branch. I don't think BigFloat should be handled separately because the upper triangular part of a LowerTriangular matrix should never be referenced. I think the solution is to change the getindex method such the indices are interchanged in zero(A.data[i,j]) such the type is fetched from the part that has the non-zeros.

mfasi · 2015-08-21T16:35:42Z

@andreasnoack Thank you, I implemented the change you suggested to getindex and both copy and scale work now. Still, every time I use the data field directly I get the nasty #undef again. I noticed that in the full! methods

function full!{T,S}(A::LowerTriangular{T,S})
    B = A.data
    tril!(B)
    B
end
function full!{T,S}(A::UnitLowerTriangular{T,S})
    B = A.data
    tril!(B)
    for i = 1:size(A,1)
        B[i,i] = 1
    end
    B
end
function full!{T,S}(A::UpperTriangular{T,S})
    B = A.data
    triu!(B)
    B
end
function full!{T,S}(A::UnitUpperTriangular{T,S})
    B = A.data
    triu!(B)
    for i = 1:size(A,1)
        B[i,i] = 1
    end
    B
end

but there might be several other places as well. It looks like we don't really want to access the data field, but I'm not sure how to overcome the issue.

andreasnoack · 2015-08-25T00:27:06Z

@mfasi It doesn't look like you have pushed your latest changes.

mfasi · 2015-08-25T06:54:22Z

@andreasnoack I hadn't because some of the tests failed. I pushed them, so that you can have a look, but the tests still fail.

mfasi · 2015-08-26T14:42:33Z

@andreasnoack I think that I might limit myself to just add copy! and scale!, which are the only methods I need. Nevertheless, I've just noticed that the issue with BigFloat causes a bug (on master):

julia> A = UpperTriangular(convert(Matrix{BigFloat},randn(5,5)));

julia> scale(A,0.5)
ERROR: BoundsError: attempt to access 5x5 Base.LinAlg.UpperTriangular{Base.MPFR.BigFloat,Array{Base.MPFR.BigFloat,2}}:
 5.203245023508947797452606209844816476106643676757812500000000000000000000000000e-01  …  9.935906547387056697573370911413803696632385253906250000000000000000000000000000e-01
 0.000000000000000000000000000000000000000000000000000000000000000000000000000000         8.351518361159303571739087601599749177694320678710937500000000000000000000000000e-01
 0.000000000000000000000000000000000000000000000000000000000000000000000000000000         1.277831052353999519155536290782038122415542602539062500000000000000000000000000    
 0.000000000000000000000000000000000000000000000000000000000000000000000000000000         1.189765274846499742977812275057658553123474121093750000000000000000000000000000    
 0.000000000000000000000000000000000000000000000000000000000000000000000000000000         5.026842940210519694588242600730154663324356079101562500000000000000000000000000e-01
  at index [2,1]
 in generic_scale! at linalg/generic.jl:20
 in scale at linalg/generic.jl:5

andreasnoack · 2015-08-26T17:52:48Z

@mfasi But copy! and scale! should also be sufficient for making copy and scale work. All take a look at the test failure.

mfasi · 2015-08-27T15:04:36Z

@andreasnoack I think that as long as we explicitly use A.data we cannot have #undefs there. Using copy built on the top of my copy! cannot work for BigFloat matrices, because the part of the matrix with #undefs is accessed by other methods (full and triu!, for istance).

My last push tries to solve the BigFloat issue by maintaining the old copy function for BigFloats and the improved version for all the other types.

andreasnoack · 2015-08-27T16:07:51Z

base/linalg/triangular.jl

+    end
+    return A
+end
+copy{S<:AbstractTriangular{BigFloat}}(A::S) = S(copy(A.data))


This shouldn't be there. If it gives an error then there is another bug to big. Special casing BigFloat is not an option. There could be many other element types with the same problem, e.g. DecFP and Matrix.

mfasi · 2015-08-28T10:19:52Z

@andreasnoack I had never thought about data types not in the test, thank you for pointing this out. I think that properly solving the issue with BigFloat (and similar data types) falls out of the scope of this PR, and I'm not sure how huge a change that could be. For the time being, I just restored the previous version of copy for triangular matrices (removing the special case for BigFloats).

As far as I can tell, this commit solves a bug (scale for triangular matrices) and adds two methods I need, without breaking anything else. I don't know whether the behaviour of #undef-using types really is a bug, but I can open an issue about it if you think it's worth it.

Add scale, scale! and copy! for triangular matrices

andreasnoack · 2015-08-28T13:20:29Z

@mfasi Thanks. If you still see visible undefs from working with Triangular matrices it would be helpful file an issue such that we have it fixed.

mfasi · 2015-08-28T13:29:44Z

@andreasnoack Thank you, I think I know a way to produce the #undef bug from a legitimate computation, I'll file an issue.

stevengj reviewed Aug 18, 2015
View reviewed changes

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

mfasi force-pushed the scale_copy_triangular branch 2 times, most recently from 7210458 to e300aa4 Compare August 20, 2015 13:14

andreasnoack reviewed Aug 21, 2015
View reviewed changes

mfasi force-pushed the scale_copy_triangular branch from e300aa4 to 7549546 Compare August 25, 2015 06:50

mfasi force-pushed the scale_copy_triangular branch 2 times, most recently from ef20861 to 3df7014 Compare August 27, 2015 13:59

andreasnoack reviewed Aug 27, 2015
View reviewed changes

Add scale, scale! and copy! for triangular matrices

4a3fa05

mfasi force-pushed the scale_copy_triangular branch from 3df7014 to 4a3fa05 Compare August 27, 2015 17:07

mfasi changed the title ~~Add scale and copy for triangular matrices~~ Add scale, scale! and copy! for triangular matrices Aug 28, 2015

andreasnoack added a commit that referenced this pull request Aug 28, 2015

Merge pull request #12683 from mfasi/scale_copy_triangular

d23449d

Add scale, scale! and copy! for triangular matrices

andreasnoack merged commit d23449d into JuliaLang:master Aug 28, 2015

mfasi deleted the scale_copy_triangular branch August 28, 2015 13:24

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Add scale, scale! and copy! for triangular matrices #12683

Add scale, scale! and copy! for triangular matrices #12683

mfasi commented Aug 18, 2015

stevengj Aug 18, 2015

mfasi commented Aug 21, 2015

andreasnoack Aug 21, 2015

mfasi Aug 21, 2015

mfasi commented Aug 21, 2015

andreasnoack commented Aug 21, 2015

mfasi commented Aug 21, 2015

andreasnoack commented Aug 25, 2015

mfasi commented Aug 25, 2015

mfasi commented Aug 26, 2015

andreasnoack commented Aug 26, 2015

mfasi commented Aug 27, 2015

andreasnoack Aug 27, 2015

mfasi commented Aug 28, 2015

andreasnoack commented Aug 28, 2015

mfasi commented Aug 28, 2015

Add scale, scale! and copy! for triangular matrices #12683

Add scale, scale! and copy! for triangular matrices #12683

Conversation

mfasi commented Aug 18, 2015

stevengj Aug 18, 2015

Choose a reason for hiding this comment

mfasi commented Aug 21, 2015

andreasnoack Aug 21, 2015

Choose a reason for hiding this comment

mfasi Aug 21, 2015

Choose a reason for hiding this comment

mfasi commented Aug 21, 2015

andreasnoack commented Aug 21, 2015

mfasi commented Aug 21, 2015

andreasnoack commented Aug 25, 2015

mfasi commented Aug 25, 2015

mfasi commented Aug 26, 2015

andreasnoack commented Aug 26, 2015

mfasi commented Aug 27, 2015

andreasnoack Aug 27, 2015

Choose a reason for hiding this comment

mfasi commented Aug 28, 2015

andreasnoack commented Aug 28, 2015

mfasi commented Aug 28, 2015