#127 fix: support cones for generic number types. #136

lrnv · 2021-06-29T14:54:07Z

Hey,

This should solve #127 , at least it allows me to run ( and obtain correct result) in the exemples from @migarstka, but also on my own application (fixing matrices that should be SDP in a certain moment problem), that could become a test.

This includes a dependency to GenericLinearAlgebra.

The code is allocating, I tried MutableArithmetics but it did not work.

This is probably suboptimal, since the ws variable is not needed anymore it could be removed to free space, but it'll make the code less readable.

Todo :

Check results on @migarstka 's and my exemple.
Eventually find a non-allocative eigen function in GenericLinearAlgebra
Find the right function in MutableArithmetics to avoid allocations.
Construct a test from my exemple & @migarstka 's one.

CLAassistant · 2021-06-29T14:54:11Z

All committers have signed the CLA.

migarstka · 2021-06-29T15:16:35Z

src/convexset.jl

@@ -203,6 +203,12 @@ function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T}
 rank_k_update!(X, ws)
 end

+function _project!(X::AbstractMatrix{BigFloat}, ws::PsdBlasWorkspace{BigFloat})
+ w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X))


This reallocates w and Z at every iteration. Is there a non-allocating version of GenericAlgebra.eigen ? eigen! perhaps?

There is a function with that name but I do not understand how it works (and where it stores the result). I tried several things and no luck yet :/

src/convexset.jl

migarstka · 2021-06-29T15:18:03Z

Thanks for submitting the PR :) I added two comments regarding the code.
I am quite busy this week, but happy to add a unit test for this next week.

lrnv · 2021-06-29T15:47:06Z

Thanks for taking a look, i changed some things. I did not manage to make GenericLinearAlgebra.eigen! work though.

blegat · 2021-06-30T11:37:36Z

src/convexset.jl

@@ -203,6 +203,11 @@ function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T}
 rank_k_update!(X, ws)
 end

+function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T}
+ w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X))
+ X .= Z'LinearAlgebra.Diagonal(max.(w,0))*Z


The matrix multiplication don't exploit the mutability of BigFloat, MutableArithmetics.jl has implementations that do

Trying MutableAritchmetics.mul!() gives ERROR: mutable_operate!(*, ::Matrix{BigFloat}, ::Diagonal{BigFloat, Vector{BigFloat}}, ::Adjoint{BigFloat, Matrix{BigFloat}}) is not implemented yet.
And MutableArithmetics.mul() just allocates as much as taking simple products.
What do you want me to do ? I'm not used to MutableArithmetics.

blegat · 2021-06-30T13:07:44Z

src/convexset.jl

@@ -205,6 +205,7 @@ end

 function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T}
 w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X))
+ #X .= MutableArithmetics.mul!(Z,LinearAlgebra.Diagonal(max.(w,0)),Z') " does not work. 


Multiplication with Diagonal is not implemented yet. You can do broadcast instead

julia> g(n) = (x = big.(ones(n)); A = big.(ones(n, n)); @time A .* x) g (generic function with 1 method) julia> j(n) = (x = big.(ones(n)); A = big.(ones(n, n)); @time MA.mul!.(A, x)) j (generic function with 1 method) julia> g(1000); 0.138854 seconds (2.00 M allocations: 114.441 MiB) julia> g(1000); 0.294522 seconds (2.00 M allocations: 114.441 MiB, 48.67% gc time) julia> g(1000); 0.139574 seconds (2.00 M allocations: 114.441 MiB) julia> j(1000); 0.182560 seconds (2 allocations: 7.629 MiB, 54.97% gc time) julia> j(1000); 0.082594 seconds (2 allocations: 7.629 MiB)

Note that here it will modify the entries of Z so of Z' as well. You can see that it still allocates 7MB because it still needs to allocate the resulting matrix.

Here is what I would do

julia> X = big.(rand(3,3)) 3×3 Matrix{BigFloat}: 0.867983 0.677627 0.457641 0.42301 0.29595 0.787097 0.531521 0.938297 0.101263 julia> MA.mutable_operate!(zero, X) julia> X 3×3 Matrix{BigFloat}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 julia> w = big.([-3.0, 0.0, 2.0]) 3-element Vector{BigFloat}: 3.0 0.0 2.0 julia> Z = BigFloat.(reshape(1:9, 3, 3)) 3×3 Matrix{BigFloat}: 1.0 4.0 7.0 2.0 5.0 8.0 3.0 6.0 9.0 julia> buffer = zero(X) 3×3 Matrix{BigFloat}: 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 julia> for i in 1:3 w[i] <= 0 && continue z = view(Z, i, :) MA.mutable_operate_to!(buffer, *, z, z') for I in eachindex(X) X[I] = MA.add_mul!(X[I], w[i], buffer[I]) end end julia> X 3×3 Matrix{BigFloat}: 21.0 48.0 75.0 48.0 120.0 192.0 75.0 192.0 309.0 julia> Z' * Diagonal(w) * Z 3×3 Matrix{BigFloat}: 21.0 48.0 75.0 48.0 120.0 192.0 75.0 192.0 309.0

Your code and outputs are not concordant, however your algorithm is OK. Sadly i needed Z*D*Z' and not Z'*D*Z, but i can translate.
Thanks :)

Looks like the code I wrote is what it should be, however i have a bug on my machine :

ERROR: MethodError: mutable_operate!(::typeof(MutableArithmetics.add_mul), ::BigFloat, ::BigFloat, ::BigFloat) is ambiguous. Candidates: mutable_operate!(op::Union{typeof(MutableArithmetics.add_mul), typeof(MutableArithmetics.sub_mul)}, x, y, z, args::Vararg{Any, N}) where N in MultivariatePolynomials at C:\Users\u009192\.julia\packages\MultivariatePolynomials\kvnmd\src\operators.jl:357 mutable_operate!(op::Function, x::BigFloat, args::Vararg{Any, N}) where N in MutableArithmetics at C:\Users\u009192\.julia\packages\MutableArithmetics\8xkW3\src\bigfloat.jl:87 Possible fix, define mutable_operate!(::Union{typeof(MutableArithmetics.add_mul), typeof(MutableArithmetics.sub_mul)}, ::BigFloat, ::Any, ::Any, ::Vararg{Any, N}) where N Stacktrace:

Looks like MultivariatePolynomials defined mutability in a way that was too greedy. @blegat what do you think ?

For other types than BigFloat, say Double64 from DoubleFloats.jl, it works correctly however.

lrnv · 2021-07-01T08:04:24Z

After @blegat fixes of type piracy at JuliaAlgebra/MultivariatePolynomials.jl#164, the code runs. The results on my use case seems to be alright for MultiFloats but not for BigFloats ! Which is crazy.

blegat · 2021-07-01T11:28:29Z

What is the problem you have for BigFloat? MultiFloats don't implement MutableArithmetics.jl

blegat · 2021-07-01T11:31:47Z

src/convexset.jl

+function _project!(X::AbstractMatrix{T}, ws::PsdBlasWorkspace{T}) where {T}
+ w,Z = GenericLinearAlgebra.eigen(GenericLinearAlgebra.Hermitian(X))
+ # The follwoing lines uses MutableArithmetics to perform the operation : X .= Z*LinearAlgebra.Diagonal(max.(w,0))*Z'
+ X = zero(X)


Because of this line, you don't modify the input X of this function

It should be MA mutable_operate!(zero, X)

Which gives :

ERROR: ArgumentError: Cannot call `mutable_operate!(zero, ::Base.ReshapedArray{BigFloat, 2, SubArray{BigFloat, 1, Vector{BigFloat}, Tuple{UnitRange{Int64}}, true}, Tuple{}})` as objects of type `Base.ReshapedArray{BigFloat, 2, SubArray{BigFloat, 1, Vector{BigFloat}, Tuple{UnitRange{Int64}}, true}, Tuple{}}` cannot be modifed to equal the result of the operation. Use `operate!` instead which returns the value of the result (possibly modifying the first argument) to write generic code that also works when the type cannot be modified.

Apparently, objects of type Base.ReshapedArray{BigFloat, 2, SubArray{BigFloat, 1, Vector{BigFloat}, Tuple{UnitRange{Int64}}, true}, Tuple{}} are not mutable, and this is the type of X.

Which means that the X that lends in the function is not mutable ?

ping @blegat, if this is still a blocker

We should implement MA.operate(zero, ...) for that type of array. IIUC, this function is expected to modify the input array. In that case, it link the local variable X to a new array and then modify it but it is not returned or anything so I don't see how that would work

Adding support for BigFLoat cones

4e135d8

removing completely the typecheck

8607081

migarstka reviewed Jun 29, 2021

View reviewed changes

src/convexset.jl Outdated Show resolved Hide resolved

lrnv added 2 commits June 29, 2021 17:30

Better dispatch

75dc3a0

Remove D

2421625

lrnv closed this Jun 29, 2021

lrnv reopened this Jun 29, 2021

blegat reviewed Jun 30, 2021

View reviewed changes

lrnv added 2 commits June 30, 2021 14:21

formula was false

337a365

Tring MutableArithmetics + small change.

3e25ee5

lrnv changed the title ~~Adding support for BigFLoat cones: solves #127~~ #127 fix: support cones for generic number types. Jun 30, 2021

blegat reviewed Jun 30, 2021

View reviewed changes

Move to mutable arithmetics for aribitrary type PSD projections

c5b1f70

lrnv mentioned this pull request Jun 30, 2021

Looks like the code I wrote is what it should be, however i have a bug on my machine : #137

Closed

blegat mentioned this pull request Jun 30, 2021

Fix type piracy with MutableArithmetics JuliaAlgebra/MultivariatePolynomials.jl#164

Merged

lrnv added 2 commits July 1, 2021 10:49

Change comment

edfcbe9

Temporary fix

ae81709

blegat reviewed Jul 1, 2021

View reviewed changes

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

#127 fix: support cones for generic number types. #136

#127 fix: support cones for generic number types. #136

lrnv commented Jun 29, 2021 •

edited

CLAassistant commented Jun 29, 2021 •

edited

migarstka Jun 29, 2021

lrnv Jun 29, 2021

migarstka commented Jun 29, 2021 •

edited

lrnv commented Jun 29, 2021 •

edited

blegat Jun 30, 2021

lrnv Jun 30, 2021

blegat Jun 30, 2021 •

edited

blegat Jun 30, 2021

lrnv Jun 30, 2021 •

edited

lrnv Jun 30, 2021

lrnv Jun 30, 2021 •

edited

lrnv commented Jul 1, 2021 •

edited

blegat commented Jul 1, 2021

blegat Jul 1, 2021

blegat Jul 1, 2021

lrnv Jul 1, 2021

lrnv Jul 1, 2021

lrnv Jul 1, 2021 •

edited

matbesancon Sep 14, 2023

blegat Sep 14, 2023

#127 fix: support cones for generic number types. #136

Are you sure you want to change the base?

#127 fix: support cones for generic number types. #136

Conversation

lrnv commented Jun 29, 2021 • edited

CLAassistant commented Jun 29, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

migarstka commented Jun 29, 2021 • edited

lrnv commented Jun 29, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

blegat Jun 30, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lrnv Jun 30, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lrnv Jun 30, 2021 • edited

Choose a reason for hiding this comment

lrnv commented Jul 1, 2021 • edited

blegat commented Jul 1, 2021

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

lrnv Jul 1, 2021 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

lrnv commented Jun 29, 2021 •

edited

CLAassistant commented Jun 29, 2021 •

edited

migarstka commented Jun 29, 2021 •

edited

lrnv commented Jun 29, 2021 •

edited

blegat Jun 30, 2021 •

edited

lrnv Jun 30, 2021 •

edited

lrnv Jun 30, 2021 •

edited

lrnv commented Jul 1, 2021 •

edited

lrnv Jul 1, 2021 •

edited