Merge pull request #3 from lsst/master

tickets/DM-19991

LICENSE.md

@@ -1,6 +1,6 @@
 MIT License
 
-Copyright (c) 2017 Peter Melchior, Fred Moolekamp
+Copyright (c) 2017-2019 Peter Melchior, Fred Moolekamp
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

README.md

@@ -26,9 +26,10 @@ We ask that any published work that utilizes this package cites:
     title="Block-simultaneous direction method of multipliers: a proximal primal-dual splitting algorithm for nonconvex problems with multiple constraints",
     journal="Optimization and Engineering",
     year="2018",
-    month="Mar",
-    day="20",
-    issn="1573-2924",
+    month="Dec",
+    volume=19,
+    issue=4,
+    pages={871-885},
     doi="10.1007/s11081-018-9380-y",
     url="https://doi.org/10.1007/s11081-018-9380-y"
     archivePrefix="arXiv",
@@ -102,15 +103,18 @@ def prox_gradf_circle(X, step):
 L = 2         # Lipschitz constant of grad f
 step_f = 1./L # maximum step size of smooth function: 1/L
 X0 = np.array([-1,0])
-X = pa.pgm(X0, prox_gradf_circle, step_f)
+# X: updated quantity
+# convergence: if iterate difference are smaller than relative error
+# error: X^{it} - X^{it-1}
+X, convergence, error = pa.pgm(X0, prox_gradf_circle, step_f)
 # or with Nesterov acceleration
-X = pa.apgm(X0, prox_gradf_circle, step_f)  
+X, convergence, error = pa.apgm(X0, prox_gradf_circle, step_f)  
 ```
 
 If the objective function is not smooth, one can use ADMM. This also allows for two functions (the objective and one constraint ) to be satisfied, but it treats them *separately*. Unlike PGM, the constraint is only met at the end of the optimization and only within some error tolerance.
 
 ```python
-X = pa.admm(X, prox_gradf, step_f, prox_circle, e_rel=1e-3, e_abs=1e-3)
+X, convergence, error = pa.admm(X, prox_gradf, step_f, prox_circle, e_rel=1e-3, e_abs=1e-3)
 ```
 
 A fully working example to demonstrate the principle of operations is [examples/parabola.py] that find the minimum of a 2D parabola under hard boundary constraints (on a shifted circle or the intersection of lines).

examples/parabola.py

@@ -3,6 +3,11 @@
 from proxmin import algorithms as pa
 from proxmin.utils import Traceback
 
+import logging
+logging.basicConfig()
+logger = logging.getLogger('proxmin')
+logger.setLevel(logging.INFO)
+
 # location of true minimum of f
 dx,dy = 1,0.5
 
@@ -132,7 +137,7 @@ def steps_f12(j=None, Xs=None):
     return slack / L
 
 
-def plotResults(tr, label, boundary=None):
+def plotResults(tr, label="", boundary=None):
     import matplotlib.pyplot as plt
     import matplotlib.patches as patches
 
@@ -148,7 +153,6 @@ def plotResults(tr, label, boundary=None):
     else:
         traj = np.dstack([tr["X",j] for j in range(tr.N)])[:,0,:]
 
-    #ax.scatter(traj[:,0], traj[:,1], s=4, c=np.arange(1,len(traj)+1), cmap='Blues_r')
     if tr.offset == 0:
         ax.plot(traj[:,0], traj[:,1], 'b.', markersize=4, ls='-')
     else:
@@ -170,7 +174,7 @@ def plotResults(tr, label, boundary=None):
     plt.show()
 
 if __name__ == "__main__":
-    xy = np.array([-1.,-1.])
+    xy0 = np.array([-1.,-1.])
     if len(sys.argv)==2:
         boundary = sys.argv[1]
         if boundary not in ["line", "circle"]:
@@ -187,56 +191,64 @@ def plotResults(tr, label, boundary=None):
 
     # PGM without boundary
     tr = Traceback()
-    x = pa.pgm(xy, prox_gradf, step_f, max_iter=max_iter, relax=1, traceback=tr)
+    xy = xy0.copy()
+    pa.pgm(xy, prox_gradf, step_f, max_iter=max_iter, relax=1, traceback=tr)
     plotResults(tr, "PGM no boundary")
 
     # PGM
     tr = Traceback()
-    x = pa.pgm(xy, prox_gradf_, step_f, max_iter=max_iter, traceback=tr)
+    xy = xy0.copy()
+    pa.pgm(xy, prox_gradf_, step_f, max_iter=max_iter, traceback=tr)
     plotResults(tr, "PGM", boundary=boundary)
 
     # APGM
     tr = Traceback()
-    x = pa.pgm(xy, prox_gradf_, step_f, max_iter=max_iter, accelerated=True,  traceback=tr)
+    xy = xy0.copy()
+    pa.pgm(xy, prox_gradf_, step_f, max_iter=max_iter, accelerated=True,  traceback=tr)
     plotResults(tr, "PGM accelerated", boundary=boundary)
 
     # ADMM
     tr = Traceback()
-    x  = pa.admm(xy, prox_gradf, step_f, prox_g, max_iter=max_iter, traceback=tr)
+    xy = xy0.copy()
+    pa.admm(xy, prox_gradf, step_f, prox_g, max_iter=max_iter, traceback=tr)
     plotResults(tr, "ADMM", boundary=boundary)
 
     # ADMM with direct constraint projection
     prox_g_direct = None
     tr = Traceback()
-    x = pa.admm(xy, prox_gradf_, step_f, prox_g_direct, max_iter=max_iter, traceback=tr)
+    xy = xy0.copy()
+    pa.admm(xy, prox_gradf_, step_f, prox_g_direct, max_iter=max_iter, traceback=tr)
     plotResults(tr, "ADMM direct", boundary=boundary)
 
     # SDMM
     M = 2
     proxs_g = [prox_g] * M # using same constraint several, i.e. M, times
     tr = Traceback()
-    x = pa.sdmm(xy, prox_gradf, step_f, proxs_g, max_iter=max_iter, traceback=tr)
+    xy = xy0.copy()
+    pa.sdmm(xy, prox_gradf, step_f, proxs_g, max_iter=max_iter, traceback=tr)
     plotResults(tr, "SDMM", boundary=boundary)
 
     # Block-SDMM
-    XY = [np.array([xy[0]]), np.array([xy[1]])]
     if boundary == "line":
         N = 2
         M1 = 7
         M2 = 2
         proxs_g = [[prox_xline]*M1, [prox_yline]*M2]
         tr = Traceback(N)
-        x = pa.bsdmm(XY, prox_gradf12, steps_f12, proxs_g, max_iter=max_iter, traceback=tr)
+        XY = [np.array([xy0[0]]), np.array([xy0[1]])]
+        pa.bsdmm(XY, prox_gradf12, steps_f12, proxs_g, max_iter=max_iter, traceback=tr)
         plotResults(tr, "bSDMM", boundary=boundary)
 
         # bSDMM with direct constraint projection
         prox_gradf12_ = partial(prox_gradf_lim12, boundary=boundary)
         prox_g_direct = None
         tr = Traceback(N)
-        x = pa.bsdmm(XY, prox_gradf12_, steps_f12, prox_g_direct, max_iter=max_iter, traceback=tr)
+        XY = [np.array([xy0[0]]), np.array([xy0[1]])]
+        pa.bsdmm(XY, prox_gradf12_, steps_f12, prox_g_direct, max_iter=max_iter, traceback=tr)
         plotResults(tr, "bSDMM direct", boundary=boundary)
 
         # BPGM
         tr = Traceback(N)
-        x = pa.bpgm(XY, prox_gradf12_, steps_f12, max_iter=max_iter, accelerated=True, traceback=tr)
+        XY = [np.array([xy0[0]]), np.array([xy0[1]])]
+        pa.bpgm(XY, prox_gradf12_, steps_f12, max_iter=max_iter, accelerated=True, traceback=tr)
         plotResults(tr, "bPGM", boundary=boundary)

examples/unmixing.py

@@ -1,12 +1,18 @@
 from proxmin import nmf
 from proxmin.utils import Traceback
+from proxmin import operators as po
 from scipy.optimize import linear_sum_assignment
 import numpy as np
 import matplotlib.pyplot as plt
 import time
-from proxmin import operators as po
 from functools import partial
 
+# initialize and run NMF
+import logging
+logging.basicConfig()
+logger = logging.getLogger('proxmin')
+logger.setLevel(logging.INFO)
+
 def generateComponent(m):
     """Creates oscillating components to be mixed"""
     freq = 25*np.random.random()
@@ -54,13 +60,12 @@ def match(A, S, trueS):
     # if noise is variable, specify variance matrix of the same shape as Y
     W = None
 
-    # initialize and run NMF
-    A0 = np.array([generateAmplitudes(k) for i in range(b)])
-    S0 = np.array([generateComponent(n) for i in range(k)])
+    A = np.array([generateAmplitudes(k) for i in range(b)])
+    S = np.array([generateComponent(n) for i in range(k)])
     p1 = partial(po.prox_unity_plus, axis=1)
     proxs_g=[[p1], None]
     tr = Traceback(2)
-    A, S = nmf(Y, A0, S0, W=W, prox_A=p1, e_rel=1e-6, e_abs=1e-6/noise**2, traceback=tr)
+    nmf(Y, A, S, W=W, prox_A=p1, e_rel=1e-6, e_abs=1e-6/noise**2, traceback=tr)
     # sort components to best match inputs
     A, S = match(A, S, trueS)