Problem solving a simple SDP in ProximalAlgorithms.jl

[Shuvomoy]

I am having issues running the following trivial problem in ProximalAlgorithms.jl:

using ProximalAlgorithms, ProximalOperators

# optimization problem: 
# minimize || X ||_{*} // nuclear norm
# subject to X >= 0

f = IndPSD()  # indicator function for PSD matrices
g = NuclearNorm() # nuclear norm function

x0 = zeros(10,10) # this is also the optimal solution

solver = ProximalAlgorithms.DouglasRachford{R}(gamma=Float64(1), tol=TOL, verbose = true, freq = 1)

y, z, it = solver(x0, f=f, g=g)

where I am getting the following error:

MethodError: no method matching prox!(::Array{Float64,2}, ::IndPSD, ::Array{Float64,2}, ::Float64)
Closest candidates are:
  prox!(::AbstractArray{T,N} where N, !Matched::CubeNormL2{R}, ::AbstractArray{T,N} where N, ::R) where {R, T<:Union{Complex{R}, R}} at C:\Users\shuvo\.julia\packages\ProximalOperators\18mAn\src\functions\cubeNormL2.jl:42
  prox!(::AbstractArray{R,N} where N, !Matched::ElasticNet{R}, ::AbstractArray{R,N} where N, ::R) where R<:Real at C:\Users\shuvo\.julia\packages\ProximalOperators\18mAn\src\functions\elasticNet.jl:39
  prox!(::AbstractArray{T,N} where N, !Matched::HuberLoss, ::AbstractArray{T,N} where N, ::Real) where T<:Union{Real, Complex} at C:\Users\shuvo\.julia\packages\ProximalOperators\18mAn\src\functions\huberLoss.jl:60
  ...
iterate(::ProximalAlgorithms.DRS_iterable{Float64,Float64,Array{Float64,2},IndPSD,NuclearNorm{Float64}}, ::ProximalAlgorithms.DRS_state{Array{Float64,2}}) at douglasrachford.jl:30
iterate(::ProximalAlgorithms.DRS_iterable{Float64,Float64,Array{Float64,2},IndPSD,NuclearNorm{Float64}}) at douglasrachford.jl:30
iterate(::ProximalAlgorithms.IterationTools.HaltingIterable{ProximalAlgorithms.DRS_iterable{Float64,Float64,Array{Float64,2},IndPSD,NuclearNorm{Float64}},ProximalAlgorithms.var"#stop#25"{ProximalAlgorithms.DouglasRachford{Float64}}}) at iterationtools.jl:20
iterate at iterators.jl:587 [inlined]
iterate at iterators.jl:585 [inlined]
iterate at iterators.jl:139 [inlined]
iterate at iterators.jl:138 [inlined]
iterate at iterationtools.jl:79 [inlined]
iterate at iterationtools.jl:79 [inlined]
iterate at iterationtools.jl:52 [inlined]
loop at iterationtools.jl:126 [inlined]
#_#24(::IndPSD, ::NuclearNorm{Float64}, ::ProximalAlgorithms.DouglasRachford{Float64}, ::Array{Float64,2}) at douglasrachford.jl:70
(::Core.var"#kw#Any")(::NamedTuple{(:f, :g),Tuple{IndPSD,NuclearNorm{Float64}}}, ::ProximalAlgorithms.DouglasRachford{Float64}, ::Array{Float64,2}) at none:0
top-level scope at testing_proximal_algorithms.jl:58

Please let me know how I can resolve the issue.

[lostella]

@Shuvomoy I think the problem is that prox! is only defined for matrices of Symmetric or Hermitian type, see its definition int he case of IndPSD. So you should be able to get your code working by doing

using LinearAlgebra

x0 = Symmetric(zeros(10, 10))

I hope that helps, let me know how that works out for you, I'm really curious!

[Shuvomoy]

Thanks so much for your quick response! Yes, I just tried your suggestion, but I am getting a new error now.

Full Code:

using ProximalAlgorithms, ProximalOperators, LinearAlgebra

# optimization problem: 
# minimize || X ||_{*} // nuclear norm
# subject to X >= 0
TOL = Real(1e-4)

T = Float64

R = real(T)

f = IndPSD()  # indicator function for PSD matrices

g = NuclearNorm() # nuclear norm function

x0 = Symmetric(zeros(10,10)) # this is also the optimal solution

solver = ProximalAlgorithms.DouglasRachford{R}(gamma=Float64(1), tol=TOL, verbose = true, freq = 1)

y, z, it = solver(x0, f=f, g=g)

Error I am getting now:

ArgumentError: Cannot set a non-diagonal index in a symmetric matrix
setindex! at symmetric.jl:225 [inlined]
_setindex! at abstractarray.jl:1104 [inlined]
setindex! at abstractarray.jl:1074 [inlined]
macro expansion at broadcast.jl:909 [inlined]
macro expansion at simdloop.jl:77 [inlined]
copyto! at broadcast.jl:908 [inlined]
copyto! at broadcast.jl:863 [inlined]
materialize! at broadcast.jl:822 [inlined]
iterate(::ProximalAlgorithms.DRS_iterable{Float64,Float64,Symmetric{Float64,Array{Float64,2}},IndPSD,NuclearNorm{Float64}}, ::ProximalAlgorithms.DRS_state{Symmetric{Float64,Array{Float64,2}}}) at douglasrachford.jl:31
iterate at douglasrachford.jl:30 [inlined]
iterate at iterationtools.jl:20 [inlined]
iterate at iterators.jl:587 [inlined]
iterate at iterators.jl:585 [inlined]
iterate at iterators.jl:139 [inlined]
iterate at iterators.jl:138 [inlined]
iterate(::ProximalAlgorithms.IterationTools.SamplingIterable{Base.Iterators.Enumerate{Base.Iterators.Take{ProximalAlgorithms.IterationTools.HaltingIterable{ProximalAlgorithms.DRS_iterable{Float64,Float64,Symmetric{Float64,Array{Float64,2}},IndPSD,NuclearNorm{Float64}},ProximalAlgorithms.var"#stop#25"{ProximalAlgorithms.DouglasRachford{Float64}}}}}}, ::Base.Iterators.Enumerate{Base.Iterators.Take{ProximalAlgorithms.IterationTools.HaltingIterable{ProximalAlgorithms.DRS_iterable{Float64,Float64,Symmetric{Float64,Array{Float64,2}},IndPSD,NuclearNorm{Float64}},ProximalAlgorithms.var"#stop#25"{ProximalAlgorithms.DouglasRachford{Float64}}}}}) at iterationtools.jl:79
iterate at iterationtools.jl:79 [inlined]
iterate at iterationtools.jl:52 [inlined]
loop(::ProximalAlgorithms.IterationTools.TeeIterable{ProximalAlgorithms.IterationTools.SamplingIterable{Base.Iterators.Enumerate{Base.Iterators.Take{ProximalAlgorithms.IterationTools.HaltingIterable{ProximalAlgorithms.DRS_iterable{Float64,Float64,Symmetric{Float64,Array{Float64,2}},IndPSD,NuclearNorm{Float64}},ProximalAlgorithms.var"#stop#25"{ProximalAlgorithms.DouglasRachford{Float64}}}}}},ProximalAlgorithms.var"#disp#26"}) at iterationtools.jl:126
#_#24 at douglasrachford.jl:70 [inlined]
(::Core.var"#kw#Any")(::NamedTuple{(:f, :g),Tuple{IndPSD,NuclearNorm{Float64}}}, ::ProximalAlgorithms.DouglasRachford{Float64}, ::Symmetric{Float64,Array{Float64,2}}) at none:0
top-level scope at douglas_rachford_modified.jl:166

Please let me know if I am missing anything.

[lostella]

I'm assuming you're using v0.3.1. The problem appears to be in this line: apparently, setting non-diagonal entries in a Symmetric object is not allowed, not even if you're actually setting the correspondent symmetric entry to the same value.

Thank you for opening this, the first solution that comes to mind is to actually fix this on the ProximalOperators side, by adding a wrapper that packs/unpacks a Matrix into a Symmetric. For now, you should be able to patch it locally in your script by adding

using ProximalOperators
using LinearAlgebra

function ProximalOperators.prox!(y::Array{R, 2}, f::IndPSD, x::Array{R, 2}, gamma) where R
    prox!(Symmetric(y), f, Symmetric(x), gamma)
end

and then directly passing matrices (not wrapped in Symmetric) as you initially did.

Let me know if that works!

[Shuvomoy]

Thanks so much for the suggestion, it works just fine now!

In case if it helps someone else, the working full code is below.

using ProximalAlgorithms, ProximalOperators, LinearAlgebra

# optimization problem:
# minimize || X ||_{*} // nuclear norm
# subject to X >= 0

TOL = Real(1e-4)

T = Float64

R = real(T)

f = IndPSD()  # indicator function for PSD matrices

g = NuclearNorm() # nuclear norm function

x0 = zeros(10,10) # this is also the optimal solution

# define a proximal operator that will work with symmetric matrix variable
function ProximalOperators.prox!(y::Array{R, 2}, f::IndPSD, x::Array{R, 2}, gamma) where R
    prox!(Symmetric(y), f, Symmetric(x), gamma)
end

solver = ProximalAlgorithms.DouglasRachford{R}(gamma=Float64(1), tol=TOL, verbose = true, freq = 1)

y, z, it = solver(x0, f=f, g=g)

Thanks again!

[lostella]

I'm glad that helped. I'll leave this open and open an issue in ProximalOperators linking to this, it makes sense to have this wrapping/unwrapping directly there. Of course there should be the complex case as well, which will use Hermitian instead of Symmetric.

[lostella]

@Shuvomoy FYI, the patch above was added to ProximalOperators, see JuliaFirstOrder/ProximalOperators.jl#101.

Closing this, thanks again for reporting the issue!

[Shuvomoy]

Great, thanks so much for your help! I am finding both ProximalAlgorithms and ProximalOperators very useful for my research. Thank you so much for creating these excellent packages!