Adds tests for the Consensus optimizer

descent/main.py

@@ -26,6 +26,19 @@ def optional_print(self, message):
             self.display.write(message + "\n")
             self.display.flush()
 
+    def add(self, operator, *args):
+        """Adds a proximal operator to the list of operators"""
+
+        if isinstance(operator, str):
+            op = getattr(proxops, operator)(*args)
+        elif isinstance(operator, proxops.ProximalOperatorBaseClass):
+            op = operator
+        else:
+            raise ValueError("operator must be a string or a subclass of ProximalOperator")
+
+        self.operators.append(op)
+        return self
+
 
 class Consensus(Optimizer):
     def __init__(self, tau=(10., 2., 2.), tol=(1e-6, 1e-3)):
@@ -47,17 +60,6 @@ def __init__(self, tau=(10., 2., 2.), tol=(1e-6, 1e-3)):
         self.tau = namedtuple('tau', ('init', 'inc', 'dec'))(*tau)
         self.tol = namedtuple('tol', ('primal', 'dual'))(*tol)
 
-    def add(self, operator, *args):
-        """Adds a proximal operator to the list of operators"""
-
-        if isinstance(operator, str):
-            op = getattr(proxops, operator)(*args)
-        elif issubclass(operator, proxops.ProximalOperatorBaseClass):
-            op = operator
-
-        self.operators.append(op)
-        return self
-
     def minimize(self, x0, display=None, maxiter=np.Inf):
 
         self.theta = x0

tests/test_proxops.py

@@ -3,41 +3,48 @@
 """
 from descent import proxops
 import numpy as np
+import pytest
 
 
 def randomseed(func):
-    """
-    Sets the seed of numpy's random number generator
-
-    (for reproducible tests)
-    """
-
+    """Sets the seed of numpy's random number generator"""
     def wrapper(*args, **kwargs):
         np.random.seed(123)
         return func(*args, **kwargs)
-
     return wrapper
 
 
+def test_baseclass():
+    with pytest.raises(TypeError):
+        op = proxops.ProximalOperatorBaseClass()
+        op(np.array([0.,]), 1.0)
+
+
 @randomseed
 def test_nucnorm():
-
     pen = 1.
     rho = 0.1
     tol = 1.
 
-    op = proxops.nucnorm(pen)
-
     X = 2 * np.outer(np.random.randn(50), np.random.randn(25))
     V = X + 0.5 * np.random.randn(50, 25)
 
-    nn = lambda A: np.linalg.svd(A, compute_uv=False).sum()
+    # compute the nuclear norm
+    def nn(A):
+        return np.linalg.svd(A, compute_uv=False).sum()
 
+    # test nucnorm
+    op = proxops.nucnorm(pen)
     assert np.abs(nn(X) - nn(op(V, rho)) - pen / rho) <= tol
 
+    # test nucnorm (w/ reshape)
+    op = proxops.nucnorm(pen, newshape=(50, 25))
+    v = op(V.ravel(), rho)
+    assert v.shape == (50*25,)
+    assert np.abs(nn(X) - nn(v.reshape(X.shape)) - pen / rho) <= tol
 
-def test_sparse():
 
+def test_sparse():
     pen = 0.1
     rho = 0.1
     v = np.linspace(-5, 5, 1e3)
@@ -51,7 +58,6 @@ def test_sparse():
 
 
 def test_nonneg():
-
     op = proxops.nonneg()
 
     v = np.array([-2., 0.54, -0.2, 24.])
@@ -63,7 +69,6 @@ def test_nonneg():
 
 @randomseed
 def test_linsys():
-
     A = np.random.randn(50,20)
     x = np.random.randn(20,)
     y = A.dot(x)
@@ -75,7 +80,6 @@ def test_linsys():
 
 
 def test_squared_error():
-
     xobs = np.array([-2., -1., 0., 1., 2.])
     tol = 1e-5
 
@@ -88,7 +92,6 @@ def test_squared_error():
 
 @randomseed
 def test_smooth():
-
     # noisy sine
     x_true = np.sin(np.linspace(0, np.pi, 100))
     x_obs = x_true + np.random.randn(x_true.size) * 0.2
@@ -118,7 +121,6 @@ def f_df(x):
 
 @randomseed
 def test_sdc():
-
     x = np.random.randn(10, 5)
     A = x.T.dot(x)
     op = proxops.sdcone()
@@ -131,3 +133,10 @@ def test_sdc():
 
     assert np.allclose(np.maximum(u0, 0), u1)
     assert np.linalg.norm(Ahat - A) <= np.linalg.norm(Aobs - A)
+
+
+def test_linear():
+    x = np.array([-2., -1., 0., 1., 2.])
+    weights = np.array([-1., 1., 0., -2., 2.])
+    op = proxops.linear(weights)
+    assert np.allclose(op(x, 0.5), x - 2 * weights)

tests/test_sparse_regression.py

@@ -2,16 +2,48 @@
 Test suite for sparse regression
 """
 import numpy as np
-from descent import algorithms
-from descent.proxops import sparse
+from descent import algorithms, Consensus
+from descent.proxops import sparse, linsys
 
 
 def generate_sparse_system(n=100, m=50, p=0.1, eta=0.05, seed=1234):
-    """
-    Generate a sparse, noisy system
+    """Generate a sparse, noisy system (Ax = y)
+
+    Parameters
+    ----------
+    n : int
+        Number of variables
+
+    m : int
+        Number of observations
+
+    p : float
+        Probability of non-zero variable
+
+    eta : float
+        Noise level (standard deviation of additive Gaussian noise)
+
+    seed : int
+        Random seed
+
+    Returns
+    -------
+    A : array_like
+        Sensing matrix (m x n)
 
+    y : array_like
+        Observations (m,)
+
+    x_true : array_like
+        True sparse signal (n,)
+
+    xls : array_like
+        Least squares solution (n,)
+
+    ls_error : float
+        Error (2-norm) of the least squares solution
     """
-    print("Generating data for sparse regression")
+
     global x_true, x_obs, A
 
     # define the seed
@@ -24,26 +56,28 @@ def generate_sparse_system(n=100, m=50, p=0.1, eta=0.05, seed=1234):
     A = np.random.randn(m, n)
     y = A.dot(x_true)
 
-    return A, y, x_true
+    # least squares solution
+    xls = np.linalg.lstsq(A, y)[0]
 
+    # least squares error
+    ls_error = np.linalg.norm(xls - x_true, 2)
 
-def test_sparse_regression():
-    """
-    Test sparse regression
+    return A, y, x_true, xls, ls_error
 
-    """
 
-    A, y, x_true = generate_sparse_system()
+def relative_error(xhat, x_true, baseline):
+    test_err = np.linalg.norm(xhat - x_true, 2)
+    return test_err / baseline
 
-    # least squares solution
-    xls = np.linalg.lstsq(A, y)[0]
+
+def test_projected_gradient_descent():
+    """Test sparse regression"""
+
+    A, y, x_true, xls, ls_error = generate_sparse_system()
 
     # helper function to test relative error
     def test_error(xhat):
-        test_err = np.linalg.norm(xhat - x_true, 2)
-        naive_err = np.linalg.norm(xls - x_true, 2)
-        err_ratio = test_err / naive_err
-        assert err_ratio <= 0.01
+        assert relative_error(xhat, x_true, ls_error) <= 0.01
 
     # Proximal gradient descent and Accelerated proximal gradient descent
     for algorithm in ['sgd', 'nag']:
@@ -57,10 +91,27 @@ def f_df(x):
 
         # optimizer
         opt = getattr(algorithms, algorithm)(lr=5e-3)
-        opt.operators.append(sparse(10.0))
+        opt.add(sparse(10.0))
 
         # run it
         res = opt.minimize(f_df, xls, display=None, maxiter=5000)
 
         # test
         test_error(res.x)
+
+
+def test_consensus():
+    """Test the consensus optimizer (ADMM)"""
+
+    A, y, x_true, xls, ls_error = generate_sparse_system()
+
+    # optimizer
+    opt = Consensus()
+    opt.add(linsys(A, y))
+    opt.add(sparse(10.0))
+
+    # run it
+    res = opt.minimize(xls, display=None, maxiter=5000)
+
+    # test
+    assert relative_error(res.x, x_true, ls_error) <= 0.05