Skip to content

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

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.

Sign up for GitHub

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

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

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

Show file tree
Hide file tree
Changes from 6 commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
4e135d8
Adding support for BigFLoat cones
lrnv Jun 29, 2021
8607081
removing completely the typecheck
lrnv Jun 29, 2021
75dc3a0
Better dispatch
lrnv Jun 29, 2021
2421625
Remove D
lrnv Jun 29, 2021
337a365
formula was false
lrnv Jun 30, 2021
3e25ee5
Tring MutableArithmetics + small change.
lrnv Jun 30, 2021
c5b1f70
Move to mutable arithmetics for aribitrary type PSD projections
lrnv Jun 30, 2021
edfcbe9
Change comment
lrnv Jul 1, 2021
ae81709
Temporary fix
lrnv Jul 1, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
3 changes: 2 additions & 1 deletion Project.toml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,6 +7,7 @@ version = "0.8.1"
AMD = "14f7f29c-3bd6-536c-9a0b-7339e30b5a3e"
COSMOAccelerators = "bbd8fffe-5ad0-4d78-a55e-85575421b4ac"
DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"
Expand All @@ -29,8 +30,8 @@ DataStructures = "^0.17.0, ^0.18.0"
IterTools = "^1"
MathOptInterface = "~0.9.16"
QDLDL = "~0.1.4"
Requires = "^1"
Reexport = "0.2, ^1"
Requires = "^1"
UnsafeArrays = "0.3, 1"
julia = "^1"

Expand Down
2 changes: 1 addition & 1 deletion src/COSMO.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
__precompile__()
module COSMO

using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg, DataStructures, Requires, Printf, IterTools
using SparseArrays, LinearAlgebra, SuiteSparse, QDLDL, Pkg, DataStructures, Requires, Printf, IterTools, GenericLinearAlgebra
using Reexport
using COSMOAccelerators

Expand Down
16 changes: 12 additions & 4 deletions src/convexset.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -183,7 +183,7 @@ for (syevr, elty) in
end #@eval
end #for

function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T}
function _project!(X::AbstractMatrix, ws::PsdBlasWorkspace{T}) where{T <: Union{Float32,Float64}}

#computes the upper triangular part of the projection of X onto the PSD cone

Expand All @@ -203,6 +203,12 @@ 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))

@migarstka migarstka Jun 29, 2021

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

@lrnv lrnv Jun 29, 2021

Copy link
Copy Markdown
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 :/

#X .= MutableArithmetics.mul!(Z,LinearAlgebra.Diagonal(max.(w,0)),Z') " does not work.

@blegat blegat Jun 30, 2021 •
edited
Loading

Copy link
Copy Markdown
Contributor

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.

@blegat blegat Jun 30, 2021

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

@lrnv lrnv Jun 30, 2021 •
edited
Loading

Copy link
Copy Markdown
Author

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

@lrnv lrnv Jun 30, 2021

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

@lrnv lrnv Jun 30, 2021 •
edited
Loading

Copy link
Copy Markdown
Author

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.

X .= Z*LinearAlgebra.Diagonal(max.(w,0))*Z'
end

function rank_k_update!(X::AbstractMatrix, ws::COSMO.PsdBlasWorkspace{T}) where {T}
n = size(X, 1)
@. X = zero(T)
Expand Down Expand Up @@ -274,9 +280,11 @@ function project!(x::AbstractVector{T}, cone::Union{PsdCone{T}, DensePsdCone{T}}
symmetrize_upper!(X)
_project!(X, cone.work)

#fill in the lower triangular part
for j=1:n, i=1:(j-1)
X[j,i] = X[i,j]
#fill in the lower triangular part, only needed when calling BLAS.
if T <: Union{Float32,Float64}
for j=1:n, i=1:(j-1)
X[j,i] = X[i,j]
end
end
end
return nothing
Expand Down
2 changes: 0 additions & 2 deletions src/interface.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -388,8 +388,6 @@ function type_checks(constraints::Vector{COSMO.Constraint{T}}) where {T <: Abstr
return nothing
end
type_checks(convex_set::AbstractConvexSet) = nothing
type_checks(convex_set::Union{PsdCone{BigFloat}, PsdConeTriangle{BigFloat}}) = throw(ArgumentError("COSMO currently does not support the combination of PSD constraints and BigFloat."))



function check_A_dim(A::Union{AbstractVector{<:Real},AbstractMatrix{<:Real}}, n::Int)
Expand Down