I am new to JuMP and trying to solve the nuclear norm minimization problem for low-rank matrix completion:
\text{min}_\mathbf{X} || \mathbf{X} ||_* \ \text{ subject to }\ x_{ij} = m_{ij},\ \ (i, j) \in \Omega
Where \Omega is the set of observed entries.
I am having trouble writing and registering the objective function in a way compatible with JuMP, and I haven’t been able to find a similar example in the tutorials or documentation. (I am using Julia 1.10.3 with JuMP v1.23.5.)
So far I have tried the following (and variations around it), but they all return a similar error:
# Create example data for the problem
using LinearAlgebra
using Random
rng = MersenneTwister(1234)
n = 20
r = 2
Y = randn(rng,n,r) * randn(rng,r,n)
X = Y .* (rand(rng, n, n) .< 1/2)
mask = X .!= 0
# Set-up the optimization problem
using JuMP, ProxSDP
model = Model(ProxSDP.Optimizer)
@variable(model, X[1:n, 1:n])
@constraint(model, X[mask] .== B[mask])
# 1st attempt:
nuclearnorm(A) = sum(LinearAlgebra.svdvals(A))
@operator(model, op_nuclearnorm, 1, nuclearnorm)
@objective(model, Min, nuclearnorm(X))
# ERROR: MethodError: no method matching svdvals!(::Matrix{NonlinearExpr})
# 2nd attempt
nuclearnorm(A) = sum(LinearAlgebra.diag(sqrt(A' * A)))
@operator(model, op_nuclearnorm, 1, nuclearnorm)
@objective(model, Min, nuclearnorm(X))
# ERROR: MethodError: no method matching sqrt(::Matrix{QuadExpr})
# 3rd attempt
f(x...) = sum(LinearAlgebra.svdvals(reshape(collect(x), n, n)))
@operator(model, myf, n*n, f)
@objective(model, Min, f(X...))
# ERROR: MethodError: no method matching svdvals!(::Matrix{NonlinearExpr})
I am guessing I need to implement a sparse version of the objective, but I am not sure how.
I have read in answers to other threads that JuMP is not really designed for structured matrix inputs; does this mean I should be using a different package altogether, or is there a way to formulate it?
(For context, in Python, the following works:
import numpy as np
import cvxpy as cp
def nuclear_norm_minimization(A, mask):
X = cp.Variable(A.shape)
objective = cp.Minimize(cp.norm(X, "nuc"))
constraints = [X[mask] == A[mask]]
problem = cp.Problem(objective, constraints)
problem.solve()
return X.value
But I would prefer to stay in Julia, and trying to see if a convenient solution exists on this side.)