Mercurial > hg > pycsalgos
changeset 3:537f7798e186
BP test working
| author | nikcleju |
|---|---|
| date | Sun, 23 Oct 2011 20:02:34 +0000 |
| parents | 735a0e24575c |
| children | 4393ad5bffc1 |
| files | pyCSalgos/BP/l1qc.py tests/BPgentest.m tests/BPtestdata.mat tests/l1qc_test.py |
| diffstat | 4 files changed, 119 insertions(+), 15 deletions(-) [+] |
line wrap: on
line diff
--- a/pyCSalgos/BP/l1qc.py Fri Oct 21 13:53:49 2011 +0000 +++ b/pyCSalgos/BP/l1qc.py Sun Oct 23 20:02:34 2011 +0000 @@ -200,7 +200,7 @@ #--------------------- # Original Matab code: - ## check if the matrix A is implicit or explicit + ## check if the mix A is implicit or explicit #largescale = isa(A,'function_handle'); # ## line search parameters @@ -330,8 +330,8 @@ fu1 = x - u fu2 = -x - u - fe = 1/2*(np.vdot(r,r) - epsilon**2) - f = u.sum() - (1/tau)*(np.log(-fu1).sum() + np.log(-fu2).sum() + math.log(-fe)) + fe = 1.0/2*(np.vdot(r,r) - epsilon**2) + f = u.sum() - (1.0/tau)*(np.log(-fu1).sum() + np.log(-fu2).sum() + math.log(-fe)) niter = 0 done = 0 @@ -363,7 +363,7 @@ #h11pfun = @(z) sigx.*z - (1/fe)*At(A(z)) + 1/fe^2*(atr'*z)*atr; h11pfun = lambda z: sigx*z - (1.0/fe)*At(A(z)) + 1.0/(fe**2)*np.dot(np.dot(atr.T,z),atr) dx,cgres,cgiter = cgsolve(h11pfun, w1p, cgtol, cgmaxiter, 0) - if (cgres > 1/2): + if (cgres > 1.0/2): print 'Cannot solve system. Returning previous iterate. (See Section 4 of notes for more information.)' xp = x.copy() up = u.copy() @@ -372,7 +372,8 @@ Adx = A(dx) else: #H11p = diag(sigx) - (1/fe)*AtA + (1/fe)^2*atr*atr'; - H11p = np.diag(sigx) - (1/fe)*AtA + (1/fe)**2*np.dot(atr,atr.T) + # Attention: atr is column vector, so atr*atr' means outer(atr,atr) + H11p = np.diag(sigx) - (1.0/fe)*AtA + (1.0/fe)**2*np.outer(atr,atr) #opts.POSDEF = true; opts.SYM = true; #[dx,hcond] = linsolve(H11p, w1p, opts); #if (hcond < 1e-14) @@ -382,7 +383,8 @@ #end try: dx = scipy.linalg.solve(H11p, w1p, sym_pos=True) - hcond = 1.0/scipy.linalg.cond(H11p) + #dx = np.linalg.solve(H11p, w1p) + hcond = 1.0/np.linalg.cond(H11p) except scipy.linalg.LinAlgError: print 'Matrix ill-conditioned. Returning previous iterate. (See Section 4 of notes for more information.)' xp = x.copy() @@ -412,9 +414,9 @@ # -fu1(ifu1)./(dx(ifu1)-du(ifu1)); -fu2(ifu2)./(-dx(ifu2)-du(ifu2)); ... # (-bqe+sqrt(bqe^2-4*aqe*cqe))/(2*aqe) # ])); - smax = min(1,np.concatenate((-fu1[ifu1]/(dx[ifu1]-du[ifu1]) , -fu2[ifu2]/(-dx[ifu2]-du[ifu2]) , (-bqe + math.sqrt(bqe**2-4*aqe*cqe))/(2*aqe)),0).min()) - - s = (0.99)*smax + smax = min(1,np.concatenate( (-fu1[ifu1]/(dx[ifu1]-du[ifu1]) , -fu2[ifu2]/(-dx[ifu2]-du[ifu2]) , np.array([ (-bqe + math.sqrt(bqe**2-4*aqe*cqe))/(2*aqe) ]) ) , 0).min()) + + s = 0.99 * smax # backtracking line search suffdec = 0 @@ -427,9 +429,9 @@ #fu1p = xp - up; fu2p = -xp - up; fep = 1/2*(rp'*rp - epsilon^2); fu1p = xp - up fu2p = -xp - up - fep = 0.5*(np.vdot(r,r) - epsilon**2) + fep = 0.5*(np.vdot(rp,rp) - epsilon**2) #fp = sum(up) - (1/tau)*(sum(log(-fu1p)) + sum(log(-fu2p)) + log(-fep)); - fp = up.sum() - (1.0/tau)*(np.log(-fu1p).sum() + np.log(-fu2p).sum() + log(-fep)) + fp = up.sum() - (1.0/tau)*(np.log(-fu1p).sum() + np.log(-fu2p).sum() + math.log(-fep)) #flin = f + alpha*s*(gradf'*[dx; du]); flin = f + alpha*s*np.dot(gradf.T , np.concatenate((dx,du),0)) #suffdec = (fp <= flin); @@ -460,7 +462,7 @@ f = fp #lambda2 = -(gradf'*[dx; du]); - lambda2 = -np.dot(gradf , np.concatenate((dx,du),0)) + lambda2 = -np.dot(gradf.T , np.concatenate((dx,du),0)) #stepsize = s*norm([dx; du]); stepsize = s * np.linalg.norm(np.concatenate((dx,du),0)) niter = niter + 1 @@ -623,7 +625,7 @@ newtonmaxiter = 50 #N = length(x0); - N = x0.size() + N = x0.size # starting point --- make sure that it is feasible if largescale: @@ -633,7 +635,7 @@ AAt = lambda z: A(At(z)) # TODO: implement cgsolve w,cgres,cgiter = cgsolve(AAt, b, cgtol, cgmaxiter, 0) - if (cgres > 1/2): + if (cgres > 1.0/2): print 'A*At is ill-conditioned: cannot find starting point' xp = x0.copy() return xp @@ -652,6 +654,7 @@ #end try: w = scipy.linalg.solve(np.dot(A,A.T), b, sym_pos=True) + #w = np.linalg.solve(np.dot(A,A.T), b) hcond = 1.0/scipy.linalg.cond(np.dot(A,A.T)) except scipy.linalg.LinAlgError: print 'A*At is ill-conditioned: cannot find starting point' @@ -673,7 +676,7 @@ # choose initial value of tau so that the duality gap after the first # step will be about the origial norm - tau = max(((2*N+1)/np.abs(x0).sum()), 1) + tau = max(((2*N+1.0)/np.abs(x0).sum()), 1) lbiter = math.ceil((math.log(2*N+1)-math.log(lbtol)-math.log(tau))/math.log(mu)) #disp(sprintf('Number of log barrier iterations = #d\n', lbiter));
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/tests/BPgentest.m Sun Oct 23 20:02:34 2011 +0000 @@ -0,0 +1,36 @@ +% Run BP and save parameters and solutions as reference test data +% to check if other algorithms are correct + +numA = 10; +numY = 100; + +sizesA{1} = [50 100]; +sizesA{2} = [20 25]; +sizesA{3} = [10 120]; +sizesA{4} = [15 100]; +sizesA{5} = [70 100]; +sizesA{6} = [80 100]; +sizesA{7} = [90 100]; +sizesA{8} = [99 100]; +sizesA{9} = [100 100]; +sizesA{10} = [250 400]; + +for i = 1:numA + sz = sizesA{i}; + cellA{i} = randn(sz); + cellY{i} = randn(sz(1), numY); + for j = 1:numY + Eps{i}(j) = rand / 5; % restrict from 0 to 20% if measurements + X0{i}(:,j) = cellA{i} \ cellY{i}(:,j); + end + %X0{i} = zeros(sz(1), numY); + +end + +for i = 1:numA + for j = 1:numY + Xr{i}(:,j) = l1qc_logbarrier_quiet(X0{i}(:,j), cellA{i}, [], cellY{i}(:,j), Eps{i}(j)); + end +end + +save BPtestdata \ No newline at end of file
--- a/tests/l1qc_test.py Fri Oct 21 13:53:49 2011 +0000 +++ b/tests/l1qc_test.py Sun Oct 23 20:02:34 2011 +0000 @@ -3,5 +3,70 @@ Created on Fri Oct 21 14:28:08 2011 @author: Nic + +Test BP algorithm """ +import numpy as np +import numpy.linalg +import scipy.io +import unittest +#import sys +#sys.path.append("D:\Nic\Dev\pyCSalgos\trunk") +#sys.path.append("D:\Nic\Dev\pyCSalgos\trunk\pyCSalgos\BP") +#import pyCSalgos +from pyCSalgos.BP.l1qc import l1qc_logbarrier +#from l1qc import l1qc_logbarrier + +class BPresults(unittest.TestCase): + def testResults(self): + mdict = scipy.io.loadmat('./data/BPtestdata.mat') + + # A = system matrix + # Y = matrix with measurements (on columns) + # X0 = matrix with initial solutions (on columns) + # Eps = vector with epsilon + # Xr = matrix with correct solutions (on columns) + for A,Y,X0,Eps,Xr in zip(mdict['cellA'].squeeze(),mdict['cellY'].squeeze(),mdict['cellX0'].squeeze(),mdict['cellEps'].squeeze(),mdict['cellXr'].squeeze()): + for i in np.arange(Y.shape[1]): + xr = l1qc_logbarrier(X0[:,i], A, np.array([]), Y[:,i], Eps.squeeze()[i]) + + # check if found solution is the same as the correct cslution + diff = numpy.linalg.norm(xr - Xr[:,i]) + err1 = numpy.linalg.norm(Y[:,i] - np.dot(A,xr)) + err2 = numpy.linalg.norm(Y[:,i] - np.dot(A,Xr[:,i])) + norm1 = numpy.linalg.norm(xr,1) + norm2 = numpy.linalg.norm(Xr[:,i],1) + print 'diff = ',diff + print 'err1 = ',err1 + print 'err2 = ',err2 + print 'norm1 = ',norm1 + print 'norm2 = ',norm2 + + # It seems Matlab's linsolve and scipy solve are slightly different + # Therefore make a more robust condition: + # OK; if solutions are close enough (diff < 1e-6) + # or + # ( + # they fulfill the constraint close enough (differr < 1e-6) + # and + # Python solution has l1 norm no more than 1e-6 larger as the reference solution + # (i.e. either norm1 < norm2 or norm1>norm2 not by more than 1e-6) + # ) + # + # ERROR: else + differr = abs((err1 - err2)) + diffnorm = norm1 - norm2 # intentionately no abs(), since norm1 < norm2 is good + if diff < 1e-6 or (differr < 1e-6 and (diffnorm < 1e-6)): + isok = True + else: + isok = False + + if not isok: + print "should raise" + #self.assertTrue(isok) + +if __name__ == "__main__": + unittest.main(verbosity=2) + #suite = unittest.TestLoader().loadTestsFromTestCase(CompareResults) + #unittest.TextTestRunner(verbosity=2).run(suite) \ No newline at end of file