Fixed nuclearnorm and sumlargesteigs for complex variables (#409)

* added tests for hcat fix

* nuclearnorm fix for complex variables

* sumlargesteigs fix for non-positive and complex variables

* Update src/problem_depot/problems/sdp.jl

Co-authored-by: Eric Hanson <[email protected]>

* Revert hcat tests which were already added

Co-authored-by: Eric Hanson <[email protected]>

* sumlargesteigs: updated comment to reflect support for non-PSD arguments

* nuclearnorm: added reference to Watrous' TQI for complex case

* added sumlargesteigs to docs/src/operations.md

* fix sumlargesteigs: sign is now NoSign

Co-authored-by: Eric Hanson <[email protected]>

Author: dstahlke

docs/src/operations.md

@@ -106,13 +106,14 @@ Semidefinite Program Representable Functions
 An optimization problem using these functions can be solved by any SDP
 solver (including SCS and Mosek).
 
-| operation          | description                       | vexity  | slope         | notes                          |
-| ------------------ | --------------------------------- | ------- | ------------- | ------------------------------ |
-| `nuclearnorm(x)`   | sum of singular values of $x$ | convex  | not monotonic | none                           |
-| `operatornorm(x)`  | max of singular values of $x$ | convex  | not monotonic | none                           |
-| `eigmax(x)`     | max eigenvalue of $x$         | convex  | not monotonic | none                           |
-| `eigmin(x)`     | min eigenvalue of $x$         | concave | not monotonic | none                           |
-| `matrixfrac(x, P)` | $x^TP^{-1}x$                  | convex  | not monotonic | IC: P is positive semidefinite |
+| operation              | description                       | vexity  | slope         | notes                          |
+| ------------------     | --------------------------------- | ------- | ------------- | ------------------------------ |
+| `nuclearnorm(x)`       | sum of singular values of $x$     | convex  | not monotonic | none                           |
+| `operatornorm(x)`      | max of singular values of $x$     | convex  | not monotonic | none                           |
+| `eigmax(x)`            | max eigenvalue of $x$             | convex  | not monotonic | none                           |
+| `eigmin(x)`            | min eigenvalue of $x$             | concave | not monotonic | none                           |
+| `matrixfrac(x, P)`     | $x^TP^{-1}x$                      | convex  | not monotonic | IC: P is positive semidefinite |
+| `sumlargesteigs(x, k)` | sum of top $k$ eigenvalues of $x$ | convex  | not monotonic | IC: P symmetric                |
 
 Exponential + SDP representable Functions
 -----------------------------------------

src/atoms/sdp_cone/nuclearnorm.jl

@@ -44,13 +44,22 @@ nuclearnorm(x::AbstractExpr) = NuclearNormAtom(x)
 #            [U A; A' V] ⪰ 0
 # see eg Recht, Fazel, Parillo 2008 "Guaranteed Minimum-Rank Solutions of Linear Matrix Equations via Nuclear Norm Minimization"
 # http://arxiv.org/pdf/0706.4138v1.pdf
+#
+# The complex case is example 1.20 of Watrous' "The Theory of Quantum Information"
+# (the operator A is negated but this doesn't affect the norm)
+# https://cs.uwaterloo.ca/~watrous/TQI/TQI.pdf
 function conic_form!(x::NuclearNormAtom, unique_conic_forms)
     if !has_conic_form(unique_conic_forms, x)
         A = x.children[1]
         m, n = size(A)
-        U = Variable(m,m)
-        V = Variable(n,n)
-        p = minimize(.5*(tr(U) + tr(V)), [U A; A' V] ⪰ 0)
+        if sign(A) == ComplexSign()
+            U = ComplexVariable(m,m)
+            V = ComplexVariable(n,n)
+        else
+            U = Variable(m,m)
+            V = Variable(n,n)
+        end
+        p = minimize(.5*real(tr(U) + tr(V)), [U A; A' V] ⪰ 0)
         cache_conic_form!(unique_conic_forms, x, p)
     end
     return get_conic_form(unique_conic_forms, x)

src/atoms/sdp_cone/sumlargesteigs.jl

@@ -1,7 +1,7 @@
 #############################################################################
 # sumlargesteigs.jl
-# Handles top k eigenvalues of a symmetric positive definite matrix
-# (and imposes the constraint that its argument be PSD)
+# Handles top k eigenvalues of a symmetric or Hermitian matrix
+# (and imposes the constraint that its argument be symmetric or Hermitian)
 # All expressions and atoms are subtypes of AbstractExpr.
 # Please read expressions.jl first.
 #############################################################################
@@ -27,7 +27,7 @@ struct SumLargestEigs <: AbstractExpr
 end
 
 function sign(x::SumLargestEigs)
-    return Positive()
+    return NoSign()
 end
 
 function monotonicity(x::SumLargestEigs)
@@ -58,11 +58,17 @@ function conic_form!(x::SumLargestEigs, unique_conic_forms)
         A = x.children[1]
         k = x.children[2]
         m, n = size(A)
-        Z = Variable(n, n)
+        if sign(A) == ComplexSign()
+            Z = ComplexVariable(n, n)
+        else
+            Z = Variable(n, n)
+        end
         s = Variable()
-        p = minimize(s*k + tr(Z),
+        # The two inequality constraints have the side effect of constraining A to be symmetric,
+        # since only symmetric matrices can be positive semidefinite.
+        p = minimize(s*k + real(tr(Z)),
                      Z + s*Matrix(1.0I, n, n) - A ⪰ 0,
-                     A ⪰ 0, Z ⪰ 0)
+                     Z ⪰ 0)
         cache_conic_form!(unique_conic_forms, x, p)
     end
     return get_conic_form(unique_conic_forms, x)

src/problem_depot/problems/sdp.jl

@@ -83,6 +83,21 @@ end
     end
 end
 
+@add_problem sdp function sdp_nuclear_norm_atom_complex(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+    A = [1 2im 3 4; 4im 3im 2 1; 4 5 6 7]
+    y = ComplexVariable(3, 4)
+    p = minimize(nuclearnorm(y), y == A; numeric_type = T)
+
+    if test
+        @test vexity(p) == ConvexVexity()
+    end
+    handle_problem!(p)
+    if test
+        @test p.optval ≈ sum(svdvals(A)) atol=atol rtol=rtol
+        @test evaluate(nuclearnorm(y)) ≈ sum(svdvals(A)) atol=atol rtol=rtol
+    end
+end
+
 @add_problem sdp function sdp_operator_norm_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable((3,3))
     p = minimize(opnorm(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12; numeric_type = T)
@@ -97,6 +112,21 @@ end
     end
 end
 
+@add_problem sdp function sdp_operator_norm_atom_complex(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
+    A = [1 2im 3 4; 4im 3im 2 1; 4 5 6 7]
+    y = ComplexVariable(3, 4)
+    p = minimize(opnorm(y), y == A; numeric_type = T)
+
+    if test
+        @test vexity(p) == ConvexVexity()
+    end
+    handle_problem!(p)
+    if test
+        @test p.optval ≈ maximum(svdvals(A)) atol=atol rtol=rtol
+        @test evaluate(opnorm(y)) ≈ maximum(svdvals(A)) atol=atol rtol=rtol
+    end
+end
+
 @add_problem sdp function sdp_sigma_max_atom(handle_problem!, ::Val{test}, atol, rtol, ::Type{T}) where {T, test}
     y = Variable((3,3))
     p = minimize(sigmamax(y), y[2,1]<=4, y[2,2]>=3, sum(y)>=12; numeric_type = T)
@@ -199,6 +229,20 @@ end
         @test p.optval ≈ 8.4853 atol=atol rtol=rtol
     end
 
+    A = [1   -2im  3   4
+         2im  7    3im 5
+         3   -3im  2   9
+         4    5    9   4]
+
+    x = ComplexVariable(4, 4)
+    p = minimize(sumlargesteigs(x, 3), x == A; numeric_type = T)
+
+    handle_problem!(p)
+
+    if test
+        @test p.optval ≈ sum(eigvals(A)[2:end]) atol=atol rtol=rtol
+    end
+
     x1 = Semidefinite(3)
     p1 = minimize(eigmax(x1), x1[1,1]>=4; numeric_type = T)