In [1]:

    

import regreg.api as rr
import numpy as np
%load_ext rmagic


    

In [2]:

    

%%R -o X,Y
library(lars)
data(diabetes)
X = diabetes$x
Y = diabetes$y


    

    






Loaded lars 1.1

In [3]:

    

X = np.hstack([X, np.ones((X.shape[0],1))])
#betah = np.dot(np.linalg.pinv(X), Y)
#Yhat = np.dot(X,betah)
#smallest_l2_bound = np.linalg.norm(Y-Yhat)
#null_soln = Y-Y.mean()
#R = smallest_l2_bound / np.linalg.norm(null_soln)


    

In [4]:

    

l2 = rr.l2norm.affine(X,-Y,bound=0.65*np.linalg.norm(Y))
l1 = rr.l1norm(X.shape[1], lagrange=1)
l2s = l2.smoothed(rr.identity_quadratic(1.e-12,0,0,0))
problem = rr.simple_problem(l2s, l1)


    

In [5]:

    

primal_nesta, dual_nesta = rr.nesta(None, l1, l2)
np.linalg.norm(Y - np.dot(X, primal_nesta)) / np.linalg.norm(Y)


    

    
Out[5]:





0.65000000000001912

In [6]:

    

l1_lagrange = rr.l1norm(X.shape[1],lagrange=np.fabs(primal_nesta).sum())
loss = rr.squared_error(X,Y, coef=2)
newsoln = rr.simple_problem(loss, l1_lagrange).solve()
np.linalg.norm(Y - np.dot(X,newsoln)) / np.linalg.norm(Y)


    

    
Out[6]:





0.31646486501863402

In [7]:

    

transform, atom = l2.dual
primal_tfocs, dual_tfocs = rr.tfocs(l1, transform, atom)
np.linalg.norm(Y - np.dot(X, primal_tfocs)) / np.linalg.norm(Y)


    

    
Out[7]:





0.65000002888829478

In [8]:

    

%timeit primal_tfocs, dual_tfocs = rr.tfocs(l1, transform, atom)


    

    




1 loops, best of 3: 2.47 s per loop

In [9]:

    

%timeit primal_nesta, dual_nesta = rr.nesta(None, l1, l2)


    

    




1 loops, best of 3: 2.24 s per loop

In [10]:

    

np.linalg.norm(primal_tfocs - primal_nesta) / (1+np.linalg.norm(primal_nesta))


    

    
Out[10]:





9.3282530767351398e-08

In [11]:

    

np.linalg.norm(dual_tfocs - dual_nesta) / (1+np.linalg.norm(dual_nesta))


    

    
Out[11]:





6.5600464362864291e-05

In [20]:

    

n, p = 200, 5000
X = np.random.standard_normal((n, p))
l1 = rr.l1norm(p, lagrange=1)
beta = np.zeros(p)
beta[:10] = 10
Y = np.dot(X, beta)

constraint = rr.zero_constraint.affine(X,-Y)
transform, atom = constraint.dual
primal_tfocs, dual_tfocs = rr.tfocs(l1, transform, atom)


    

In [21]:

    

np.linalg.norm(Y - np.dot(X, primal_tfocs)) / np.linalg.norm(Y)


    

    
Out[21]:





1.7792174654304693e-09

In [21]:

    




    

In [22]:

    

primal_tfocs[:20]


    

    
Out[22]:





array([ 10.        ,  10.        ,  10.00000001,  10.        ,
         9.99999998,   9.99999999,  10.00000001,   9.99999997,
        10.        ,  10.00000002,   0.        ,  -0.        ,
         0.        ,   0.        ,   0.        ,   0.        ,
         0.        ,  -0.        ,  -0.        ,   0.        ])

In [ ]: