some test problems in sandbox

Author: dance858

cvxpy/sandbox/tests_problem_references/NMF_fro.py

@@ -0,0 +1,16 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import cvxpy as cp
+np.random.seed(0)
+
+n, m, k = 5, 4, 2
+noise_level = 0.05
+X_true = np.random.rand(n, k)
+Y_true = np.random.rand(k, m)
+A_noisy = X_true @ Y_true + noise_level * np.random.randn(n, m)
+A_noisy = np.clip(A_noisy, 0, None)
+X = cp.Variable((n, k), nonneg=True)
+Y = cp.Variable((k, m), nonneg=True)
+obj = cp.sum_squares(A_noisy - X @ Y)
+problem = cp.Problem(cp.Minimize(obj))
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)

cvxpy/sandbox/tests_problem_references/boolean_ls.py

@@ -0,0 +1,15 @@
+import numpy as np 
+import cvxpy as cp 
+import pdb
+
+n, m = 100, 25
+A = np.random.randn(m, n)
+b = np.random.randn(m)
+
+x = cp.Variable(n)
+obj = cp.Minimize(cp.sum_squares(A @ x - b))
+constraints = [cp.square(x) == 4]
+problem = cp.Problem(obj, constraints)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+pdb.set_trace()

cvxpy/sandbox/tests_problem_references/phase_ret_cmplex.py

@@ -0,0 +1,21 @@
+import numpy as np
+import cvxpy as cp
+import pdb
+np.random.seed(0)
+
+# Problem dimensions
+m, n = 50, 10
+x_true = 10 * (np.random.randn(n) + 1j*np.random.randn(n))
+A = np.random.randn(m, n) + 1j*np.random.randn(m, n)
+y = np.abs(A @ x_true)
+
+x = cp.Variable(n, complex=True)
+x0 = np.random.rand(n) + 1j*np.random.rand(n)
+x.value = x0 
+obj = cp.sum_squares(cp.abs(A @ x) - y)
+problem = cp.Problem(cp.Minimize(obj))
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+print("x0: ", x0)
+print("true x: ", x_true)
+print("recovered x: ", x.value)

cvxpy/sandbox/tests_problem_references/phase_ret_nonsmooth_real.py

@@ -0,0 +1,32 @@
+import numpy as np
+import cvxpy as cp
+import pdb
+np.random.seed(0)
+
+m, n = 300, 50
+x_true = np.random.randn(n)
+A = np.random.randn(m, n)
+y = np.abs(A @ x_true) 
+x = cp.Variable(n)
+
+t = cp.Variable(m)
+v = cp.Variable(m)
+
+#obj = cp.norm1(cp.square(A @ x) - y ** 2)
+obj = cp.sum(t)
+constraints = [cp.square(t) == cp.square(v), v == cp.square(A @ x) - y ** 2, t >= 0]
+problem = cp.Problem(cp.Minimize(obj), constraints)
+
+# With this initialization IPOPT terminates after one iteration
+x.value = 0 * np.ones(n)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_star_one = x.value
+
+# With this initialization IPOPT converges to the x_true
+x.value = np.ones(n)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_star_two = x.value
+
+print("true x:                                 ", np.round(x_true, 2))
+print("recovered x with zero initialization:   ", np.round(x_star_one, 2))
+print("recovered x with small initialization:  ", np.round(x_star_two, 2))

cvxpy/sandbox/tests_problem_references/phase_ret_real_data.py

@@ -0,0 +1,31 @@
+import numpy as np
+import cvxpy as cp
+import pdb
+np.random.seed(0)
+
+# Problem dimensions
+m, n = 50, 10
+x_true = 10 * np.random.randn(n)
+A = np.random.randn(m, n)
+y = np.abs(A @ x_true) #+  np.random.randn(m)
+
+x = cp.Variable(n)
+# if I initialize from 0 IPOPT terminates after one iteration...
+# this is extremely weird! must look into this
+x0 =  2 * x_true
+x.value = x0 
+#obj = cp.sum_squares(cp.square(A @ x) - y ** 2)
+#obj = cp.sum_squares(A @ x - y ** 2)
+#pdb.set_trace()
+#obj = cp.sum_squares(cp.abs(A @ x) - y)
+t = cp.Variable(m)
+v = cp.Variable(m)
+#constraints = [t == cp.abs(v), v == A @ x]
+constraints = [cp.square(t) == cp.square(v), v == A @ x, t >= 0]
+obj = cp.sum_squares(t - y) 
+problem = cp.Problem(cp.Minimize(obj), constraints)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+print("x0: ", x0)
+print("true x: ", x_true)
+print("recovered x: ", x.value)

cvxpy/sandbox/tests_problem_references/phase_retrieval_real.py

@@ -0,0 +1,28 @@
+import numpy as np
+import cvxpy as cp
+import pdb
+np.random.seed(0)
+
+m, n = 10, 5
+x_true = np.random.randn(n)
+A = np.random.randn(m, n)
+y = np.abs(A @ x_true) 
+x = cp.Variable(n)
+
+
+obj = cp.sum_squares(cp.square(A @ x) - y ** 2)
+problem = cp.Problem(cp.Minimize(obj))
+
+# With this initilaization IPOPT terminates after one iteration
+x.value = 0 * np.ones(n)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_star_one = x.value
+
+# With this initialization IPOPT converges to the x_true
+x.value = 1e-8 * np.ones(n)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_star_two = x.value
+
+print("true x:                                 ", np.round(x_true, 2))
+print("recovered x with zero initialization:   ", np.round(x_star_one, 2))
+print("recovered x with small initialization:  ", np.round(x_star_two, 2))

cvxpy/sandbox/tests_problem_references/powell1.py

@@ -0,0 +1,56 @@
+import numpy as np
+import pandas as pd
+import pdb 
+
+import cvxpy as cp
+
+
+x = cp.Variable(4, name='x')
+F = cp.hstack([x[0] + 10 * x[1],
+               np.sqrt(5) * (x[2] - x[3]),
+               cp.square(x[1] - 2 * x[2]),
+               np.sqrt(10) * (x[0] - x[3])])
+
+
+# formulation 1 - works (95 iterations)
+objective = cp.sum_squares(F)
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([3, -1, 0, 1])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+# formulation 2 - works (41 iterations)
+objective = cp.square(x[0] + 10 * x[1]) + \
+            cp.square(np.sqrt(5) * (x[2] - x[3])) + \
+            cp.power(x[1] - 2 * x[2], 4) + \
+            cp.square(np.sqrt(10) * (x[0] - x[3]))
+
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([3, -1, 0, 1])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+
+# formulation 3 - works (18 iterations)
+objective = cp.square(x[0] + 10 * x[1]) + \
+            cp.square(np.sqrt(5) * (x[2] - x[3])) + \
+            cp.square(cp.square(x[1] - 2 * x[2])) + \
+            cp.square(np.sqrt(10) * (x[0] - x[3]))
+
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([3, -1, 0, 1])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+# formulation 4 - works (18 iterations)
+objective = cp.power(x[0] + 10 * x[1], 2) + \
+            cp.power(np.sqrt(5) * (x[2] - x[3]), 2) + \
+            cp.square(cp.square(x[1] - 2 * x[2])) + \
+            cp.power(np.sqrt(10) * (x[0] - x[3]), 2)
+
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([3, -1, 0, 1])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+# formulation 5 - canonicalization fails
+#objective = cp.norm(F)
+#problem = cp.Problem(cp.Minimize(objective))
+#problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+#pdb.set_trace()

cvxpy/sandbox/tests_problem_references/powell2.py

@@ -0,0 +1,24 @@
+import numpy as np
+import pandas as pd
+import pdb 
+
+import cvxpy as cp
+
+
+x = cp.Variable(2, name='x')
+F = cp.hstack([1e4*x[0]*x[1] - 1,
+               cp.exp(-x[0]) + cp.exp(-x[1]) - 1.001])
+
+# formulation 1 
+objective = cp.sum_squares(F)
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([0, 1])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+# formulation 2 ()
+objective = cp.square(1e4*x[0]*x[1] - 1) + cp.square(cp.exp(-x[0]) + cp.exp(-x[1]) - 1.001)
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([0, 1])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+
+pdb.set_trace()

cvxpy/sandbox/tests_problem_references/power_system_estimation.py

@@ -0,0 +1,36 @@
+import numpy as np 
+import cvxpy as cp 
+import pdb
+
+m = 100
+n = 5
+Ps = []
+Vs = []
+x_true = np.random.randn(n)
+for i in range(m):
+    P = np.random.randn(n, n)
+    P = P + P.T
+    Ps.append(P)
+    v = x_true.T @ P @ x_true
+    Vs.append(v)
+
+x = cp.Variable(n)
+obj = 0
+for i in range(m):
+    obj += cp.square(cp.quad_form(x, Ps[i]) - Vs[i])
+
+problem = cp.Problem(cp.Minimize(obj))
+
+# With this initilaization IPOPT terminates after one iteration
+x.value = 0 * np.ones(n)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_star_one = x.value
+
+# With this initialization IPOPT converges to the x_true
+x.value = 1e-8 * np.ones(n)
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_star_two = x.value
+
+print("true x:                                 ", np.round(x_true, 2))
+print("recovered x with zero initialization:   ", np.round(x_star_one, 2))
+print("recovered x with small initialization:  ", np.round(x_star_two, 2))

cvxpy/sandbox/tests_problem_references/roth.py

@@ -0,0 +1,25 @@
+import numpy as np
+import pandas as pd
+import cvxpy as cp
+
+x = cp.Variable(2, name='x')
+F = cp.hstack([-13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1],
+               -29 + x[0] + ((x[1] + 1) * x[1] - 14) * x[1]])
+
+# formulation 1 
+objective = cp.sum_squares(F)
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([0.5, -2])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_opt_one = x.value
+
+# formulation 2
+objective = cp.square(-13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1]) + \
+            cp.square(-29 + x[0] + ((x[1] + 1) * x[1] - 14) * x[1])
+problem = cp.Problem(cp.Minimize(objective))
+x.value = np.array([0.5, -2])
+problem.solve(solver=cp.IPOPT, nlp=True, verbose=True)
+x_opt_two = x.value
+
+print(x_opt_one)
+print(x_opt_two)