Define copy!, setindex!, and arithmetic for Symmetric and Hermitian #19228

ararslan · 2016-11-05T03:43:22Z

This implements copy! for Symmetric and Hermitian types, as well as + and - for Symmetric and Hermitian with UniformScaling.

tkelman · 2016-11-05T14:20:53Z

base/linalg/symmetric.jl

@@ -133,6 +133,8 @@ convert{T}(::Type{AbstractMatrix{T}}, A::Hermitian) = Hermitian(convert(Abstract

 copy{T,S}(A::Symmetric{T,S}) = (B = copy(A.data); Symmetric{T,typeof(B)}(B,A.uplo))
 copy{T,S}(A::Hermitian{T,S}) = (B = copy(A.data); Hermitian{T,typeof(B)}(B,A.uplo))
+copy!(dest::Symmetric, src::Symmetric) = copy!(dest.data, src.data)
+copy!(dest::Hermitian, src::Hermitian) = copy!(dest.data, src.data)


this is not correct if uplo differs

Good catch! Will fix, thanks!

should write a test too that would have caught that

tkelman · 2016-11-05T14:23:47Z

base/linalg/uniformscaling.jl

@@ -56,6 +56,10 @@ ishermitian(J::UniformScaling) = isreal(J.λ)
 (-)(B::BitArray{2}, J::UniformScaling)      = Array(B) - J
 (-)(J::UniformScaling, B::BitArray{2})      = J - Array(B)

+for t in (:Symmetric, :Hermitian), op in (:+, :-)
+    @eval ($op)(X::$t, J::UniformScaling) = ($t)(($op)(X.data, J), Symbol(X.uplo))


should test that Hermitian plus a complex UniformScaling with nonzero imag component currently should error

Done. Thanks!

ararslan · 2016-11-05T19:27:20Z

While I'm at it, would it make sense to implement other operators on Symmetric/Hermitian with UniformScaling, e.g. * and /? If so, which operators?

stevengj · 2016-11-05T21:00:10Z

* and / are not nearly so useful, because you can just multiply or divide by a scalar instead.

ararslan · 2016-11-05T21:04:26Z

It appears that operations on Symmetric and Hermitian with scalars similarly fail. I could define those as well...

stevengj · 2016-11-05T23:07:29Z

That sounds useful.

ararslan · 2016-11-05T23:12:52Z

Okay, can do.

andreasnoack · 2016-11-05T12:21:53Z

base/linalg/symmetric.jl

@@ -133,6 +133,8 @@ convert{T}(::Type{AbstractMatrix{T}}, A::Hermitian) = Hermitian(convert(Abstract

 copy{T,S}(A::Symmetric{T,S}) = (B = copy(A.data); Symmetric{T,typeof(B)}(B,A.uplo))
 copy{T,S}(A::Hermitian{T,S}) = (B = copy(A.data); Hermitian{T,typeof(B)}(B,A.uplo))
+copy!(dest::Symmetric, src::Symmetric) = copy!(dest.data, src.data)


The copy method should return the wrapped version, i.e. something like (copy!(dest.data, src.data);dest)

andreasnoack · 2016-11-05T12:22:50Z

base/linalg/uniformscaling.jl

@@ -56,6 +56,10 @@ ishermitian(J::UniformScaling) = isreal(J.λ)
 (-)(B::BitArray{2}, J::UniformScaling)      = Array(B) - J
 (-)(J::UniformScaling, B::BitArray{2})      = J - Array(B)

+for t in (:Symmetric, :Hermitian), op in (:+, :-)
+    @eval ($op)(X::$t, J::UniformScaling) = ($t)(($op)(X.data, J), Symbol(X.uplo))


Did you try to define a setindex! method for Symmetric first? I think it might be sufficient.

Yep, that seems to have done the trick. I was able to remove all of the special methods for Symmetric and Hermitian with UniformScaling. Thanks!

tkelman · 2016-11-06T01:18:31Z

Symmetric is a view wrapper type though, the inactive triangle should never be read from so its contents don't need to be symmetric

andreasnoack · 2016-11-06T23:23:01Z

base/linalg/symmetric.jl

@@ -107,6 +107,24 @@ end
    r
 end

+function setindex!(A::Symmetric, v, i::Integer, j::Integer)
+    if (A.uplo == 'U' && i >= j) || (A.uplo == 'L' && i <= j)


I don't think this is right. It should only be allowed to modify the diagonal. Modifying off-diagonals will effectively change two values at a time so I think it should be an error.

That makes sense. Then I'm back to defining specific methods for things, since only being able to set the diagonal doesn't allow more general operations, but maybe that's fine.

ararslan · 2016-11-07T00:03:10Z

Okay, I've improved my initial implementation of copy!, restricted setindex! to the diagonals only, defined +, -, *, and / for Symmetric and Hermitian with scalars, and added more tests. The tests pass locally, so hopefully CI likes it too.

Thanks so much for your help, everyone!

andreasnoack · 2016-11-07T00:08:16Z

base/linalg/symmetric.jl

@@ -218,6 +248,10 @@ A_mul_B!{T<:BlasComplex,S<:StridedMatrix}(C::StridedMatrix{T}, A::StridedMatrix{
 *(A::HermOrSym, B::HermOrSym) = full(A)*full(B)
 *(A::StridedMatrix, B::HermOrSym) = A*full(B)

+for T in (:Symmetric, :Hermitian), op in (:+, :-, :*, :/)
+    @eval ($op)(A::$T, x::Number) = ($T)(($op)(A.data, x), Symbol(A.uplo))


The Hermitian operations should probably be restricted to operations with x::Real.

I had thought about that, but there are cases where operating with a complex value could still produce a real diagonal. But restricting to Real does seem the safer option.

tkelman · 2016-11-07T00:50:16Z

base/linalg/symmetric.jl

+
+function setindex!(A::Hermitian, v, i::Integer, j::Integer)
+    i == j || throw(ArgumentError("Cannot set a non-diagonal index in a Hermitian matrix"))
+    isreal(v) || throw(ArgumentError("Cannot set a nonreal value to the diagonal in a Hermitian matrix"))


little backwards? cannot set a diagonal entry to a non-real value?

Good call. Fixed. Thanks!

ararslan · 2016-11-07T03:12:24Z

The CI failure appears to be unrelated. (One of the tests in the ambiguity test set is returning a non-boolean value.)

tkelman · 2016-11-07T03:24:42Z

That's a messy hack so when that test fails it tells you what the ambiguities were. But it's related.

ararslan · 2016-11-07T03:52:40Z

Ohhhh okay, so the Tuple{Method,Method} thing in the log is the list of ambiguities. I guess it's angry about my + and - definitions being ambiguous with AbstractArray{Bool,N<:Any}. Hmm.

ararslan · 2016-11-07T04:15:59Z

Any idea what I could do to ameliorate that? As long as Bool <: Number the solution is not clear to me.

tkelman · 2016-11-07T04:30:47Z

Probably have to define some ambiguity-resolving methods for the Bool case where the existing method and yours intersect?

ararslan · 2016-11-07T05:47:47Z

That seems to have done the trick. Thanks so much, @tkelman!

ararslan · 2016-11-07T18:15:28Z

Looking better now, @andreasnoack?

andreasnoack · 2016-11-08T00:25:49Z

test/linalg/symmetric.jl

@@ -250,3 +250,36 @@ let a = randn(2,2)
    cc = copy(c)
    @test conj!(c) == conj(Array(c))
 end
+
+# 19225
+let X = sparse([1 -1; -1 1])


I think these should to into the sparse tests for modularity. However, you might want to test the functionality with dense arrays as well.

I removed the sparse wrapper for these tests, then copied a simplified version of what's currently here to the sparse tests. Let me know if the format looks okay.

andreasnoack

See the comment at line 255

ararslan · 2016-11-08T23:37:20Z

All good now?

ararslan · 2016-11-10T18:43:31Z

Bump. @andreasnoack

…uliaLang#19228) * Define copy! and setindex! for Symmetric and Hermitian * Define arithmetic operations on Symmetric/Hermitian with Number * Add a test for the ambiguous case * Separate sparse and dense tests

KristofferC · 2017-04-07T18:17:33Z

base/linalg/symmetric.jl

@@ -107,6 +107,21 @@ end
    r
 end

+function setindex!(A::Symmetric, v, i::Integer, j::Integer)
+    i == j || throw(ArgumentError("Cannot set a non-diagonal index in a symmetric matrix"))


This change doesn't make sense to me. Why shouldn't you be able to set an off diagonal element?

Then it wouldn't be symmetric anymore.

This is a symmetric wrapper so it is implicitly symmetric no matter what you do with it. Setting A[i,j] = v is the same as A[i,j] = A[j,i] = v.

The problem then is that the behavior of setindex! differs from everywhere else in that it sets multiple values at once.

Yeah. We discussed that and it was considered too surprising. Are you considering a generic case where modifying the underlying buffer would be inconvenient?

I'm not sure what you mean by the right and wrong sides of a diagonal. Doesn't that only apply to triangular matrices?

I mean the .uplo field.

Only being able to assign to the upper or lower triangle depending on the value of uplo also breaks the abstraction; the user shouldn't have to care about which triangle is stored since all they want to know about it is that it's symmetric.

we could introduce a symset! function, or something similar, that generically changes values in a way that maintains symmetry.

The usefulness of a symset! not withstanding, we could also think about permitting setindex! operations that maintain symmetry. E.g. A[:,:] = 0 might be expected to work. Also A[[CartesianIndex(1,2), CartesianIndex(2,1)]] = 1, although symset! would of course be nicer for this case.

fredrikekre · 2017-06-22T06:34:56Z

Just leaving a note here that this broke the following:

julia> A = rand(2, 2); S = Symmetric(A); D = Diagonal(rand(2));

julia> S*D
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
 [1] setindex!(::Symmetric{Float64,Array{Float64,2}}, ::Float64, ::Int64, ::Int64) at ./linalg/symmetric.jl:126
 [2] scale!(::Symmetric{Float64,Array{Float64,2}}, ::Symmetric{Float64,Array{Float64,2}}, ::Array{Float64,1}) at ./linalg/matmul.jl:19
 [3] *(::Symmetric{Float64,Array{Float64,2}}, ::Diagonal{Float64}) at ./linalg/diagonal.jl:152

julia> D*S
ERROR: ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
Stacktrace:
 [1] setindex!(::Symmetric{Float64,Array{Float64,2}}, ::Float64, ::Int64, ::Int64) at ./linalg/symmetric.jl:126
 [2] scale!(::Symmetric{Float64,Array{Float64,2}}, ::Array{Float64,1}, ::Symmetric{Float64,Array{Float64,2}}) at ./linalg/matmul.jl:34
 [3] *(::Diagonal{Float64}, ::Symmetric{Float64,Array{Float64,2}}) at ./linalg/diagonal.jl:154

which I stumbled upon trying to add tests for #22396

tkelman reviewed Nov 5, 2016

View reviewed changes

ararslan force-pushed the aa/copy-symmetric branch from eb2b3d7 to 1029227 Compare November 5, 2016 17:54

andreasnoack requested changes Nov 5, 2016

View reviewed changes

ararslan force-pushed the aa/copy-symmetric branch from 1029227 to 93a1f58 Compare November 6, 2016 22:45

andreasnoack reviewed Nov 6, 2016

View reviewed changes

ararslan force-pushed the aa/copy-symmetric branch from cc19d10 to 864aa15 Compare November 7, 2016 00:01

andreasnoack reviewed Nov 7, 2016

View reviewed changes

tkelman reviewed Nov 7, 2016

View reviewed changes

Define copy! and setindex! for Symmetric and Hermitian

374cc7b

ararslan force-pushed the aa/copy-symmetric branch from 864aa15 to 832416c Compare November 7, 2016 02:29

ararslan added 2 commits November 6, 2016 20:47

Define arithmetic operations on Symmetric/Hermitian with Number

91d28bd

Add a test for the ambiguous case

3009d49

ararslan force-pushed the aa/copy-symmetric branch from 832416c to 3009d49 Compare November 7, 2016 04:54

ararslan changed the title ~~Define copy! for Symmetric and Hermitian~~ Define copy!, setindex!, and arithmetic for Symmetric and Hermitian Nov 7, 2016

andreasnoack reviewed Nov 8, 2016

View reviewed changes

andreasnoack requested changes Nov 8, 2016

View reviewed changes

Separate sparse and dense tests

de7d086

andreasnoack approved these changes Nov 10, 2016

View reviewed changes

andreasnoack merged commit 1c12e9e into JuliaLang:master Nov 10, 2016

ararslan deleted the aa/copy-symmetric branch November 10, 2016 18:47

ararslan mentioned this pull request Mar 4, 2017

fix broadcast! where the destination is a structured matrix #20886

Merged

KristofferC reviewed Apr 7, 2017

View reviewed changes

This was referenced Jun 22, 2017

carry out (c)transpose in Ax_mul_Bx methods for Hermitian and Symmetric #22396

Merged

fix matrix multiplication interaction with HermOrSym and Diagonal #22474

Merged

dalum mentioned this pull request Oct 3, 2018

Preserve types when adding/subtracting Herm/Sym/UniformScaling #29500

Merged

fredrikekre mentioned this pull request Aug 26, 2019

Allow changing arbitrary entries in symmetric/Hermitian matrices #33071

Closed

fredrikekre mentioned this pull request Oct 16, 2020

Symmertic/Hermitian Matrixes are too fussy about mutating operations JuliaLang/LinearAlgebra.jl#773

Open

Timeroot mentioned this pull request Jan 15, 2022

Allow modifying off-diagonals of SymTridiagonal #43818

Closed

jishnub mentioned this pull request Jun 9, 2023

Make SymTridiagonal off-diagonals mutable #50114

Closed

LebedevRI mentioned this pull request Jun 25, 2024

Set inequality syntax suggestion jump-dev/JuMP.jl#3770

Closed

Define copy!, setindex!, and arithmetic for Symmetric and Hermitian #19228

Define copy!, setindex!, and arithmetic for Symmetric and Hermitian #19228

Conversation

ararslan commented Nov 5, 2016

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

ararslan commented Nov 5, 2016

stevengj commented Nov 5, 2016

ararslan commented Nov 5, 2016 • edited Loading

stevengj commented Nov 5, 2016

ararslan commented Nov 5, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tkelman commented Nov 6, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

ararslan commented Nov 7, 2016

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

ararslan commented Nov 7, 2016

tkelman commented Nov 7, 2016

ararslan commented Nov 7, 2016 • edited Loading

ararslan commented Nov 7, 2016

tkelman commented Nov 7, 2016

ararslan commented Nov 7, 2016

ararslan commented Nov 7, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

andreasnoack left a comment

Choose a reason for hiding this comment

ararslan commented Nov 8, 2016

ararslan commented Nov 10, 2016

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KristofferC Apr 7, 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

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

fredrikekre commented Jun 22, 2017

ararslan commented Nov 5, 2016 •

edited

Loading

ararslan commented Nov 7, 2016 •

edited

Loading

KristofferC Apr 7, 2017 •

edited

Loading