In [3]:

    

import numpy as np
from cvxopt import solvers, lapack, matrix, spmatrix


    

In [7]:

    

solvers.options['show_progress'] = 1

m = 6 # Number of data points
n = 1 # Number of variables
tt = np.arange(1.0, m+1)
yy = np.array([8.49, 20.05, 50.65, 72.19, 129.85, 171.56])
yy[3] += 100
varphi = np.power(tt, 2)/2
varphi.shape = (m, 1)

A = matrix(varphi)
b = matrix(yy)

# Least squares solution.
xls = np.dot( np.linalg.inv(np.dot(varphi.T, varphi)), np.dot(varphi.T, yy) )
# Or
xls = np.linalg.lstsq(varphi, yy)
# Or
xls = +b
lapack.gels(+A, xls)
xls = xls[:1]

# Robust least squares.
#
# minimize  sum( h( A*x-b ))
#
# where h(u) = u^2           if |u| <= 1.0
#            = 2*(|u| - 1.0) if |u| > 1.0.
#
# Solve as a QP (see exercise 4.5):
#
# minimize    (1/2) * u'*u + 1'*v
# subject to  -u - v <= A*x-b <= u + v
#             0 <= u <= 1
#             v >= 0
#
# Variables  x (n), u (m), v(m)

novars = n+2*m
P = spmatrix([],[],[], (novars, novars))
P[n:n+m,n:n+m] = spmatrix(1.0, range(m), range(m))
q = matrix(0.0, (novars,1))
q[-m:] = 1.0

G = spmatrix([], [], [], (5*m, novars))
h = matrix(0.0, (5*m,1))

# A*x - b <= u+v
G[:m,:n] = A
G[:m,n:n+m] = spmatrix(-1.0, range(m), range(m))
G[:m,n+m:] = spmatrix(-1.0, range(m), range(m))
h[:m] = b

# -u - v <= A*x - b
G[m:2*m,:n] = -A
G[m:2*m,n:n+m] = spmatrix(-1.0, range(m), range(m))
G[m:2*m,n+m:] = spmatrix(-1.0, range(m), range(m))
h[m:2*m] = -b

# u >= 0
G[2*m:3*m,n:n+m] = spmatrix(-1.0, range(m), range(m))

# u <= 1
G[3*m:4*m,n:n+m] = spmatrix(1.0, range(m), range(m))
h[3*m:4*m] = 1.0

# v >= 0
G[4*m:,n+m:] = spmatrix(-1.0, range(m), range(m))

xh = solvers.qp(P, q, G, h)['x'][:n]

print("Huber solution: g_H = %f" % xh[0])
print("Least squares solution: g_LS = %f" % xls[0])


    

    




     pcost       dcost       gap    pres   dres
 0: -3.1983e+00  1.5683e+04  2e+05  1e+00  3e+02
 1:  4.9162e+02 -1.5174e+03  3e+03  1e-02  3e+00
 2:  3.1218e+02 -1.2353e+02  4e+02  2e-04  5e-02
 3:  1.2092e+02  7.4295e+01  5e+01  2e-05  5e-03
 4:  1.1501e+02  1.0120e+02  1e+01  4e-06  1e-03
 5:  1.1047e+02  1.0855e+02  2e+00  5e-07  1e-04
 6:  1.0996e+02  1.0968e+02  3e-01  4e-09  1e-06
 7:  1.0985e+02  1.0984e+02  2e-02  9e-11  2e-08
 8:  1.0984e+02  1.0984e+02  2e-03  6e-14  1e-11
 9:  1.0984e+02  1.0984e+02  2e-04  1e-16  3e-12
10:  1.0984e+02  1.0984e+02  3e-05  1e-16  2e-12
Optimal solution found.
Huber solution: g_H = 10.347738
Least squares solution: g_LS = 11.184167

In [9]:

    

import pylab
pylab.figure(1,facecolor='w')
pylab.plot(tt,yy,'o',
        tt, varphi*xh[0], '-g',
        tt, varphi*xls[0], '--r',
        markerfacecolor='w', markeredgecolor='b')
#pylab.axis([-11, 11, -20, 25])
pylab.xlabel('t')
pylab.ylabel('f(t)')
pylab.title('Robust regression')
pylab.show()


    

In [32]:

    

a = +b


    

In [33]:

    

a


    

    
Out[33]:





<6x1 matrix, tc='d'>

In [34]:

    

a[2]


    

    
Out[34]:





50.65

In [35]:

    

a[2] += 3


    

In [37]:

    

b[2]


    

    
Out[37]:





50.65

In [ ]: