|
|
'''
|
|
|
Created on May 30, 2014
|
|
|
|
|
|
@author: Yolian Amaro
|
|
|
'''
|
|
|
|
|
|
from irls_dn import *
|
|
|
from scipy.optimize import brentq
|
|
|
|
|
|
def irls_dn2(A,b,p,G):
|
|
|
|
|
|
# Minimize ||u||_p subject to ||A*u-b||_2^2 <= G (0 < p <= 1)
|
|
|
|
|
|
# What this function actually does is finds the lambda1 so that the solution
|
|
|
# to the following problem satisfies ||A*u-b||_2^2 <= G:
|
|
|
# Minimize lambda1*||u||_p + ||A*u-b||_2
|
|
|
|
|
|
# Start with a large lambda1, and do a line search until fidelity <= G.
|
|
|
# (Inversions with large lambda1 are really fast anyway).
|
|
|
|
|
|
# Then spin up fsolve to localize the root even better
|
|
|
|
|
|
# Line Search
|
|
|
|
|
|
alpha = 2.0; # Line search parameter
|
|
|
lambda1 = 1e5; # What's a reasonable but safe initial guess?
|
|
|
u = irls_dn(A,b,p,lambda1);
|
|
|
|
|
|
|
|
|
fid = norm(np.dot(A,u)-b)**2;
|
|
|
|
|
|
print '----------------------------------\n';
|
|
|
|
|
|
while (fid >= G):
|
|
|
lambda1 = lambda1 / alpha; # Balance between speed and accuracy
|
|
|
u = irls_dn(A,b,p,lambda1);
|
|
|
fid = norm(np.dot(A,u)-b)**2;
|
|
|
print 'lambda = %2e \t' % lambda1, '||A*u-b||^2 = %.1f' % fid;
|
|
|
|
|
|
# Refinement using brentq (equiv to fzero in matlab)
|
|
|
lambda0 = np.array([lambda1,lambda1*alpha]); # interval with zero-crossing
|
|
|
|
|
|
def myfun(lambda1):
|
|
|
return norm(np.dot(A, irls_dn(A,b,p,lambda1)) - b)**2 - G;
|
|
|
|
|
|
# Find zero-crossing at given interval (lambda1, lambda1*alpha)
|
|
|
lambda1 = brentq(myfun, lambda0[0], lambda0[1], xtol=0.01*lambda1)
|
|
|
|
|
|
u = irls_dn(A,b,p,lambda1);
|
|
|
|
|
|
return u;
|
|
|
|
|
|
|
