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

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

Open
wants to merge 9 commits into
base: master
Choose a base branch
from

Conversation

lrnv
Copy link

@lrnv lrnv commented Jun 29, 2021

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
Copy link

CLAassistant commented Jun 29, 2021

CLA assistant check
All committers have signed the CLA.

@@ -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))
Copy link
Member

Choose a reason for hiding this comment

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

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

Copy link
Author

Choose a reason for hiding this comment

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

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 Outdated Show resolved Hide resolved
@migarstka
Copy link
Member

migarstka commented Jun 29, 2021

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
Copy link
Author

lrnv commented Jun 29, 2021

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

@lrnv lrnv closed this Jun 29, 2021
@lrnv lrnv reopened this Jun 29, 2021
src/convexset.jl Outdated
@@ -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
Copy link
Contributor

Choose a reason for hiding this comment

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

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

Copy link
Author

Choose a reason for hiding this comment

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

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.

@lrnv lrnv changed the title Adding support for BigFLoat cones: solves #127 #127 fix: support cones for generic number types. Jun 30, 2021
src/convexset.jl Outdated
@@ -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.
Copy link
Contributor

@blegat blegat Jun 30, 2021

Choose a reason for hiding this comment

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

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.

Copy link
Contributor

Choose a reason for hiding this comment

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

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

Copy link
Author

@lrnv lrnv Jun 30, 2021

Choose a reason for hiding this comment

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

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 :)

Copy link
Author

Choose a reason for hiding this comment

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

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 ?

Copy link
Author

@lrnv lrnv Jun 30, 2021

Choose a reason for hiding this comment

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

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

@lrnv
Copy link
Author

lrnv commented Jul 1, 2021

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
Copy link
Contributor

blegat commented Jul 1, 2021

What is the problem you have for BigFloat? MultiFloats don't implement MutableArithmetics.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)
Copy link
Contributor

Choose a reason for hiding this comment

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

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

Copy link
Contributor

Choose a reason for hiding this comment

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

It should be MA mutable_operate!(zero, X)

Copy link
Author

Choose a reason for hiding this comment

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

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.

Copy link
Author

Choose a reason for hiding this comment

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

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.

Copy link
Author

@lrnv lrnv Jul 1, 2021

Choose a reason for hiding this comment

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

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

Copy link
Contributor

Choose a reason for hiding this comment

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

ping @blegat, if this is still a blocker

Copy link
Contributor

Choose a reason for hiding this comment

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

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

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.

5 participants