Implement bound constraints (e.g. with ADMM or Dog Leg) #3
Description
Implementing bound constraints with ADMM should be relatively easy:
http://stanford.edu/~eryu/nnlsqr.html
https://web.stanford.edu/~boyd/papers/pdf/admm_slides.pdf (Slide 27)
Quick and dirty python implementation:
import numpy as np
import scipy.optimize
def lsmr_with_init(A,b,x0):
r0 = b - scipy.sparse.linalg.aslinearoperator(A).matvec(x0)
deltax_pack = scipy.sparse.linalg.lsmr(A,r0)
return x0 + deltax_pack[0]
def admm_nnlsqr(A,b):
# ADMM parameter
# Optimal choice product of the maximum and minimum singular values
# Heuristic choice: mean of singular values, i.e. squared Frobenius norm over dimension
rho = np.dot(A.ravel(),A.ravel()) / A.shape[1]
sqrt_half_rho = np.sqrt(rho/2)
print 'sqrt_half_rho=',sqrt_half_rho
# initialisation
zero_init = False
if zero_init:
x_pack = (np.zeros(A.shape[1]),0)
z = np.zeros(A.shape[1])
u = np.zeros(A.shape[1])
else:
# from x=z=u=0 the first iteration is a normally a projection of
# the unconstrained damped least squares. Let's forget the damping and check
# whether we need to constrain things
x_pack = scipy.sparse.linalg.lsmr(A,b)
if np.all(x_pack[0]>0):
# no constraint needed
return x_pack[0]
z = np.clip(x_pack[0],0,np.inf)
u = x_pack[0]-z
# construct the damped least squares matrix
# todo try and use directly the damped version of lsmr
Adamped = scipy.sparse.linalg.LinearOperator(
A.shape+np.array([A.shape[1], 0]),
matvec = lambda y: np.concatenate([ A.dot(y), sqrt_half_rho*y ]),
rmatvec = lambda y: y[0:A.shape[0]].dot(A) + sqrt_half_rho*y[A.shape[0]:],
)
tol = 1e-5
max_iter = 5000
diff = np.inf
iter = 0
while ( diff>tol and iter<max_iter ):
iter += 1
xold = x_pack[0]
#x_pack = scipy.sparse.linalg.lsmr( Adamped, np.concatenate([ b, sqrt_half_rho*(z-u) ]) )
x_pack = (lsmr_with_init( Adamped, np.concatenate([ b, sqrt_half_rho*(z-u) ]), xold ), 0)
zold = z
z = np.clip(x_pack[0]+u,0,np.inf)
#diff = np.linalg.norm(z-zold)
#diff = np.linalg.norm(x_pack[0]-xold)
diff = np.linalg.norm(x_pack[0]-z)
u += x_pack[0]-z
return z