update of jupyter book
This commit is contained in:
Morten Hjorth-Jensen
@@ -408,7 +408,7 @@ Using the SVD we can obtain the pseudoinverse of a matrix $\bm{A}$ (labeled here
|
||||
\bm{A}_{\mathrm{PI}}= \bm{V}\bm{D}_{\mathrm{PI}}\bm{U}^T,
|
||||
\]
|
||||
!et
|
||||
where $\bm{D}_{\mathrm{PI}}$ can be calculated by creating a diagonal matrix from $\bm{Sigma}$ where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.
|
||||
where $\bm{D}_{\mathrm{PI}}$ can be calculated by creating a diagonal matrix from $\bm{\Sigma}$ where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.
|
||||
|
||||
|
||||
!bc pycod
|
||||
@@ -562,14 +562,6 @@ and using the orthogonality of the matrix $\bm{U}$ we have
|
||||
!et
|
||||
We define $\bm{\Sigma}^T\bm{\Sigma}=\tilde{\bm{\Sigma}}^2$ which is a diagonal matrix containing only the singular values squared. It has dimensionality $p \times p$.
|
||||
|
||||
This means, using the orthogonality of $\bm{V}$, that we get
|
||||
|
||||
!bt
|
||||
\[
|
||||
\bm{X}^T\bm{X}=\tilde{\bm{\Sigma}}^2.
|
||||
\]
|
||||
!et
|
||||
|
||||
We can now insert the result for the matrix $\bm{X}^T\bm{X}$ into our equation for ordinary least squares where
|
||||
|
||||
!bt
|
||||
@@ -581,10 +573,10 @@ and using our SVD decomposition of $\bm{X}$ we have
|
||||
|
||||
!bt
|
||||
\[
|
||||
\tilde{y}_{\mathrm{OLS}}=\bm{U}\bm{\Sigma}\bm{V}^T\tilde{\bm{\Sigma}}^{-2}\bm{V}\bm{\Sigma}^T\bm{U}^T\bm{y},
|
||||
\tilde{y}_{\mathrm{OLS}}=\bm{U}\bm{\Sigma}\bm{V}^T\left(\bm{V}\tilde{\bm{\Sigma}}^{2}(\bm{V}^T\right)^{-1}\bm{V}\bm{\Sigma}^T\bm{U}^T\bm{y},
|
||||
\]
|
||||
!et
|
||||
which gives us, using the orthogonality of the matrices $\bm{U}$ and $\bm{V}$,
|
||||
which gives us, using the orthogonality of the matrices $\bm{U}$ and $\bm{V}$,,
|
||||
|
||||
!bt
|
||||
\[
|
||||
@@ -592,6 +584,10 @@ which gives us, using the orthogonality of the matrices $\bm{U}$ and $\bm{V}$,
|
||||
\]
|
||||
!et
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
Note here that when we perform the multiplication of the various matrices, the orthogonal vectors of the matrix $\bm{U}$
|
||||
!bt
|
||||
\[
|
||||
@@ -1014,7 +1010,10 @@ which is just
|
||||
\end{bmatrix},
|
||||
\]
|
||||
!et
|
||||
where we wrote $$\bm{C}[\bm{x}_0,\bm{x}_1] = \bm{C}[\bm{x}]$$ to indicate that this is the covariance of the vectors $\bm{x}$ of the design/feature matrix $\bm{X}$.
|
||||
|
||||
where we wrote $\bm{C}[\bm{x}_0,\bm{x}_1]=\bm{C}[\bm{x}]$ to indicate
|
||||
that this is the covariance of the vectors $\bm{x}$ of the
|
||||
design/feature matrix $\bm{X}$.
|
||||
|
||||
It is easy to generalize this to a matrix $\bm{X}\in {\mathbb{R}}^{n\times p}$.
|
||||
|
||||
@@ -1782,282 +1781,6 @@ series of lectures.
|
||||
|
||||
|
||||
As a small addendum, we note that you can also solve this problem using the convex optimization package "CVXOPT":"https://cvxopt.org/examples/mlbook/l1regls.html". This requires, in addition to having installed _CVXOPT_, you need to download the file *l1regl.py*.
|
||||
The following code example solves the simpler problem we discussed above, where we have added the latter python file.
|
||||
|
||||
!bc pycod
|
||||
from cvxopt import matrix, spdiag, mul, div, sqrt, normal, setseed
|
||||
from cvxopt import blas, lapack, solvers, sparse, spmatrix
|
||||
import math
|
||||
|
||||
try:
|
||||
import mosek
|
||||
import sys
|
||||
__MOSEK = True
|
||||
except: __MOSEK = False
|
||||
|
||||
if __MOSEK:
|
||||
|
||||
def l1regls_mosek(A, b):
|
||||
"""
|
||||
|
||||
Returns the solution of l1-norm regularized least-squares problem
|
||||
|
||||
minimize || A*x - b ||_2^2 + e'*u
|
||||
|
||||
subject to -u <= x <= u
|
||||
|
||||
"""
|
||||
|
||||
m, n = A.size
|
||||
|
||||
env = mosek.Env()
|
||||
task = env.Task(0,0)
|
||||
task.set_Stream(mosek.streamtype.log, lambda x: sys.stdout.write(x))
|
||||
|
||||
task.appendvars( 2*n) # number of variables
|
||||
task.appendcons( 2*n) # number of constraints
|
||||
|
||||
# input quadratic objective
|
||||
Q = matrix(0.0, (n,n))
|
||||
blas.syrk(A, Q, alpha = 2.0, trans='T')
|
||||
|
||||
I = []
|
||||
for i in range(n):
|
||||
I.extend(range(i,n))
|
||||
|
||||
J = []
|
||||
for i in range(n):
|
||||
J.extend((n-i)*[i])
|
||||
|
||||
task.putqobj(I, J, list(Q[matrix(I) + matrix(J)*n]))
|
||||
task.putclist(range(2*n), list(-2*A.T*b) + n*[1.0]) # setup linear objective
|
||||
|
||||
# input constraint matrix row by row
|
||||
for i in range(n):
|
||||
task.putarow( i, [i, n+i], [1.0, -1.0])
|
||||
task.putarow( n+i, [i, n+i], [1.0, 1.0])
|
||||
|
||||
# setup bounds on constraints
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
0, n, n*[mosek.boundkey.up], n*[0.0], n*[0.0])
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
n, 2*n, n*[mosek.boundkey.lo], n*[0.0], n*[0.0])
|
||||
|
||||
# setup variable bounds
|
||||
task.putboundslice(mosek.accmode.var,
|
||||
0, 2*n, 2*n*[mosek.boundkey.fr], 2*n*[0.0], 2*n*[0.0])
|
||||
|
||||
# optimize the task
|
||||
task.putobjsense(mosek.objsense.minimize)
|
||||
task.optimize()
|
||||
task.solutionsummary(mosek.streamtype.log)
|
||||
x = n*[0.0]
|
||||
task.getsolutionslice(mosek.soltype.itr, mosek.solitem.xx, 0, n, x)
|
||||
|
||||
return matrix(x)
|
||||
|
||||
def l1regls_mosek2(A, b):
|
||||
"""
|
||||
|
||||
Returns the solution of l1-norm regularized least-squares problem
|
||||
|
||||
minimize w'*w + e'*u
|
||||
|
||||
subject to -u <= x <= u
|
||||
|
||||
A*x - w = b
|
||||
|
||||
"""
|
||||
|
||||
m, n = A.size
|
||||
|
||||
env = mosek.Env()
|
||||
task = env.Task(0,0)
|
||||
task.set_Stream(mosek.streamtype.log, lambda x: sys.stdout.write(x))
|
||||
|
||||
task.appendvars(2*n + m) # number of variables
|
||||
task.appendcons(2*n + m) # number of constraints
|
||||
|
||||
# input quadratic objective
|
||||
task.putqobj(range(2*n,2*n+m), range(2*n,2*n+m), m*[2.0])
|
||||
|
||||
task.putclist(range(2*n+m), n*[0.0] + n*[1.0] + m*[0.0]) # setup linear objective
|
||||
|
||||
# input constraint matrix row by row
|
||||
for i in range(n):
|
||||
task.putarow( i, [i, n+i], [1.0, -1.0])
|
||||
task.putarow( n+i, [i, n+i], [1.0, 1.0])
|
||||
|
||||
for i in range(m):
|
||||
task.putarow( 2*n+i, range(n) + [2*n+i], list(A[i,:]) + [-1.0])
|
||||
|
||||
# setup bounds on constraints
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
0, n, n*[mosek.boundkey.up], n*[0.0], n*[0.0])
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
n, 2*n, n*[mosek.boundkey.lo], n*[0.0], n*[0.0])
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
2*n, 2*n+m, m*[mosek.boundkey.fx], list(b), list(b))
|
||||
|
||||
# setup variable bounds
|
||||
task.putboundslice(mosek.accmode.var, 0, 2*n+m, (2*n+m)*[mosek.boundkey.fr],
|
||||
(2*n+m)*[0.0], (2*n+m)*[0.0])
|
||||
|
||||
# optimize the task
|
||||
task.putobjsense(mosek.objsense.minimize)
|
||||
task.optimize()
|
||||
task.solutionsummary(mosek.streamtype.log)
|
||||
x = n*[0.0]
|
||||
task.getsolutionslice(mosek.soltype.itr, mosek.solitem.xx, 0, n, x)
|
||||
|
||||
return matrix(x)
|
||||
|
||||
def l1regls(A, b):
|
||||
"""
|
||||
|
||||
Returns the solution of l1-norm regularized least-squares problem
|
||||
|
||||
minimize || A*x - b ||_2^2 + || x ||_1.
|
||||
|
||||
"""
|
||||
|
||||
m, n = A.size
|
||||
q = matrix(1.0, (2*n,1))
|
||||
q[:n] = -2.0 * A.T * b
|
||||
|
||||
def P(u, v, alpha = 1.0, beta = 0.0 ):
|
||||
"""
|
||||
v := alpha * 2.0 * [ A'*A, 0; 0, 0 ] * u + beta * v
|
||||
"""
|
||||
v *= beta
|
||||
v[:n] += alpha * 2.0 * A.T * (A * u[:n])
|
||||
|
||||
|
||||
def G(u, v, alpha=1.0, beta=0.0, trans='N'):
|
||||
"""
|
||||
v := alpha*[I, -I; -I, -I] * u + beta * v (trans = 'N' or 'T')
|
||||
"""
|
||||
|
||||
v *= beta
|
||||
v[:n] += alpha*(u[:n] - u[n:])
|
||||
v[n:] += alpha*(-u[:n] - u[n:])
|
||||
|
||||
h = matrix(0.0, (2*n,1))
|
||||
|
||||
|
||||
# Customized solver for the KKT system
|
||||
#
|
||||
# [ 2.0*A'*A 0 I -I ] [x[:n] ] [bx[:n] ]
|
||||
# [ 0 0 -I -I ] [x[n:] ] = [bx[n:] ].
|
||||
# [ I -I -D1^-1 0 ] [zl[:n]] [bzl[:n]]
|
||||
# [ -I -I 0 -D2^-1 ] [zl[n:]] [bzl[n:]]
|
||||
#
|
||||
# where D1 = W['di'][:n]**2, D2 = W['di'][:n]**2.
|
||||
#
|
||||
# We first eliminate zl and x[n:]:
|
||||
#
|
||||
# ( 2*A'*A + 4*D1*D2*(D1+D2)^-1 ) * x[:n] =
|
||||
# bx[:n] - (D2-D1)*(D1+D2)^-1 * bx[n:] +
|
||||
# D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] -
|
||||
# D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:]
|
||||
#
|
||||
# x[n:] = (D1+D2)^-1 * ( bx[n:] - D1*bzl[:n] - D2*bzl[n:] )
|
||||
# - (D2-D1)*(D1+D2)^-1 * x[:n]
|
||||
#
|
||||
# zl[:n] = D1 * ( x[:n] - x[n:] - bzl[:n] )
|
||||
# zl[n:] = D2 * (-x[:n] - x[n:] - bzl[n:] ).
|
||||
#
|
||||
# The first equation has the form
|
||||
#
|
||||
# (A'*A + D)*x[:n] = rhs
|
||||
#
|
||||
# and is equivalent to
|
||||
#
|
||||
# [ D A' ] [ x:n] ] = [ rhs ]
|
||||
# [ A -I ] [ v ] [ 0 ].
|
||||
#
|
||||
# It can be solved as
|
||||
#
|
||||
# ( A*D^-1*A' + I ) * v = A * D^-1 * rhs
|
||||
# x[:n] = D^-1 * ( rhs - A'*v ).
|
||||
|
||||
S = matrix(0.0, (m,m))
|
||||
Asc = matrix(0.0, (m,n))
|
||||
v = matrix(0.0, (m,1))
|
||||
|
||||
def Fkkt(W):
|
||||
|
||||
# Factor
|
||||
#
|
||||
# S = A*D^-1*A' + I
|
||||
#
|
||||
# where D = 2*D1*D2*(D1+D2)^-1, D1 = d[:n]**-2, D2 = d[n:]**-2.
|
||||
|
||||
d1, d2 = W['di'][:n]**2, W['di'][n:]**2
|
||||
|
||||
# ds is square root of diagonal of D
|
||||
ds = math.sqrt(2.0) * div( mul( W['di'][:n], W['di'][n:]),
|
||||
sqrt(d1+d2) )
|
||||
d3 = div(d2 - d1, d1 + d2)
|
||||
|
||||
# Asc = A*diag(d)^-1/2
|
||||
Asc = A * spdiag(ds**-1)
|
||||
|
||||
# S = I + A * D^-1 * A'
|
||||
blas.syrk(Asc, S)
|
||||
S[::m+1] += 1.0
|
||||
lapack.potrf(S)
|
||||
|
||||
def g(x, y, z):
|
||||
|
||||
x[:n] = 0.5 * ( x[:n] - mul(d3, x[n:]) +
|
||||
mul(d1, z[:n] + mul(d3, z[:n])) - mul(d2, z[n:] -
|
||||
mul(d3, z[n:])) )
|
||||
x[:n] = div( x[:n], ds)
|
||||
|
||||
# Solve
|
||||
#
|
||||
# S * v = 0.5 * A * D^-1 * ( bx[:n] -
|
||||
# (D2-D1)*(D1+D2)^-1 * bx[n:] +
|
||||
# D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] -
|
||||
# D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:] )
|
||||
|
||||
blas.gemv(Asc, x, v)
|
||||
lapack.potrs(S, v)
|
||||
|
||||
# x[:n] = D^-1 * ( rhs - A'*v ).
|
||||
blas.gemv(Asc, v, x, alpha=-1.0, beta=1.0, trans='T')
|
||||
x[:n] = div(x[:n], ds)
|
||||
|
||||
# x[n:] = (D1+D2)^-1 * ( bx[n:] - D1*bzl[:n] - D2*bzl[n:] )
|
||||
# - (D2-D1)*(D1+D2)^-1 * x[:n]
|
||||
x[n:] = div( x[n:] - mul(d1, z[:n]) - mul(d2, z[n:]), d1+d2 )\
|
||||
- mul( d3, x[:n] )
|
||||
|
||||
# zl[:n] = D1^1/2 * ( x[:n] - x[n:] - bzl[:n] )
|
||||
# zl[n:] = D2^1/2 * ( -x[:n] - x[n:] - bzl[n:] ).
|
||||
z[:n] = mul( W['di'][:n], x[:n] - x[n:] - z[:n] )
|
||||
z[n:] = mul( W['di'][n:], -x[:n] - x[n:] - z[n:] )
|
||||
|
||||
return g
|
||||
|
||||
return solvers.coneqp(P, q, G, h, kktsolver = Fkkt)['x'][:n]
|
||||
|
||||
|
||||
!ec
|
||||
|
||||
|
||||
Then we call the above functions and solve the problem, as done here
|
||||
|
||||
!bc pycod
|
||||
from cvxopt import matrix, normal
|
||||
|
||||
X = matrix( [ [ 2, 0, 1], [0, 1, 3]])
|
||||
y = matrix( [4, 2, 3])
|
||||
x = l1regls(X,y)
|
||||
!ec
|
||||
|
||||
_More text will be added to this example._
|
||||
|
||||
===== Linking the regression analysis with a statistical interpretation =====
|
||||
|
||||
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
After Width: | Height: | Size: 9.1 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 11 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 14 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 10 KiB |
File diff suppressed because it is too large
Load Diff
@@ -478,7 +478,8 @@ const thebe_selector_output = ".output, .cell_output"
|
||||
|
||||
<div>
|
||||
|
||||
<div class="tex2jax_ignore mathjax_ignore section" id="ridge-and-lasso-regression">
|
||||
<!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
|
||||
doconce format html chapter2.do.txt --><div class="tex2jax_ignore mathjax_ignore section" id="ridge-and-lasso-regression">
|
||||
<h1><span class="section-number">4. </span>Ridge and Lasso Regression<a class="headerlink" href="#ridge-and-lasso-regression" title="Permalink to this headline">¶</a></h1>
|
||||
<div class="section" id="mathematical-interpretation-of-ordinary-least-squares">
|
||||
<h2><span class="section-number">4.1. </span>Mathematical Interpretation of Ordinary Least Squares<a class="headerlink" href="#mathematical-interpretation-of-ordinary-least-squares" title="Permalink to this headline">¶</a></h2>
|
||||
@@ -850,7 +851,7 @@ It is used for the calculation of the inverse for singular or near singular matr
|
||||
\[
|
||||
\boldsymbol{A}_{\mathrm{PI}}= \boldsymbol{V}\boldsymbol{D}_{\mathrm{PI}}\boldsymbol{U}^T,
|
||||
\]</div>
|
||||
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{D}_{\mathrm{PI}}\)</span> can be calculated by creating a diagonal matrix from <span class="math notranslate nohighlight">\(\boldsymbol{Sigma}\)</span> where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.</p>
|
||||
<p>where <span class="math notranslate nohighlight">\(\boldsymbol{D}_{\mathrm{PI}}\)</span> can be calculated by creating a diagonal matrix from <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}\)</span> where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.</p>
|
||||
<div class="cell docutils container">
|
||||
<div class="cell_input docutils container">
|
||||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">import</span> <span class="nn">numpy</span> <span class="k">as</span> <span class="nn">np</span>
|
||||
@@ -983,11 +984,6 @@ decomposition of the design matrix.</p>
|
||||
\boldsymbol{X}^T\boldsymbol{X}=\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}\boldsymbol{V}^T.
|
||||
\]</div>
|
||||
<p>We define <span class="math notranslate nohighlight">\(\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}=\tilde{\boldsymbol{\Sigma}}^2\)</span> which is a diagonal matrix containing only the singular values squared. It has dimensionality <span class="math notranslate nohighlight">\(p \times p\)</span>.</p>
|
||||
<p>This means, using the orthogonality of <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>, that we get</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\boldsymbol{X}^T\boldsymbol{X}=\tilde{\boldsymbol{\Sigma}}^2.
|
||||
\]</div>
|
||||
<p>We can now insert the result for the matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}^T\boldsymbol{X}\)</span> into our equation for ordinary least squares where</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
@@ -996,9 +992,9 @@ decomposition of the design matrix.</p>
|
||||
<p>and using our SVD decomposition of <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span> we have</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T\tilde{\boldsymbol{\Sigma}}^{-2}\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{y},
|
||||
\tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T\left(\boldsymbol{V}\tilde{\boldsymbol{\Sigma}}^{2}(\boldsymbol{V}^T\right)^{-1}\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{y},
|
||||
\]</div>
|
||||
<p>which gives us, using the orthogonality of the matrices <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> and <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>,</p>
|
||||
<p>which gives us, using the orthogonality of the matrices <span class="math notranslate nohighlight">\(\boldsymbol{U}\)</span> and <span class="math notranslate nohighlight">\(\boldsymbol{V}\)</span>,,</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
\tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{U}^T\boldsymbol{y}=\sum_{i=0}^{p-1}\boldsymbol{u}_i\boldsymbol{u}^T_j\boldsymbol{y},
|
||||
@@ -1214,10 +1210,10 @@ covariance matrix through the <strong>np.linalg.eig()</strong> function.</p>
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>-0.13876586436927824
|
||||
3.722047011333792
|
||||
[[ 1.233528 3.58428804]
|
||||
[ 3.58428804 11.47942814]]
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>-0.01591355407242949
|
||||
3.6808538439837775
|
||||
[[ 0.96390357 2.99157584]
|
||||
[ 2.99157584 10.31120247]]
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1254,10 +1250,10 @@ a more brute force way. Here we scale the mean values for each column of the des
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.08464758160254343
|
||||
1.8503720991789538
|
||||
[[1. 0.65626043]
|
||||
[0.65626043 1. ]]
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>0.08612280083325631
|
||||
1.6149274949460215
|
||||
[[1. 0.66934291]
|
||||
[0.66934291 1. ]]
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1287,30 +1283,30 @@ this matrix we easily see that it is a positive definite matrix.</p>
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[-1.41876267 -4.93248252]
|
||||
[ 1.83687444 5.28097861]
|
||||
[ 0.37429133 0.59766 ]
|
||||
[ 0.59159438 1.71869727]
|
||||
[-0.80315282 -0.89348922]
|
||||
[-0.38748219 -2.12288563]
|
||||
[-2.08917679 -5.64933923]
|
||||
[ 0.27803645 0.89944994]
|
||||
[ 1.23703839 3.2321528 ]
|
||||
[ 0.38073947 1.86925797]]
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[[-1.04649105 -2.92658312]
|
||||
[ 0.45985488 1.43876695]
|
||||
[-0.41081513 -1.96426825]
|
||||
[ 1.75703965 4.88736621]
|
||||
[ 1.02698605 3.59304008]
|
||||
[-0.71348713 -1.98249059]
|
||||
[-0.22685646 0.37866422]
|
||||
[-0.90559087 -1.31597731]
|
||||
[-0.60349429 -3.47245463]
|
||||
[ 0.66285434 1.36393643]]
|
||||
0 1
|
||||
0 -1.418763 -4.932483
|
||||
1 1.836874 5.280979
|
||||
2 0.374291 0.597660
|
||||
3 0.591594 1.718697
|
||||
4 -0.803153 -0.893489
|
||||
5 -0.387482 -2.122886
|
||||
6 -2.089177 -5.649339
|
||||
7 0.278036 0.899450
|
||||
8 1.237038 3.232153
|
||||
9 0.380739 1.869258
|
||||
0 -1.046491 -2.926583
|
||||
1 0.459855 1.438767
|
||||
2 -0.410815 -1.964268
|
||||
3 1.757040 4.887366
|
||||
4 1.026986 3.593040
|
||||
5 -0.713487 -1.982491
|
||||
6 -0.226856 0.378664
|
||||
7 -0.905591 -1.315977
|
||||
8 -0.603494 -3.472455
|
||||
9 0.662854 1.363936
|
||||
0 1
|
||||
0 1.000000 0.977418
|
||||
1 0.977418 1.000000
|
||||
0 1.000000 0.948641
|
||||
1 0.948641 1.000000
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1367,37 +1363,37 @@ this matrix we easily see that it is a positive definite matrix.</p>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span> 0 1 2 3 4 5 6 7 \
|
||||
0 0.0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
|
||||
1 0.0 0.081253 0.083015 0.080297 0.079157 0.078040 0.071112 0.069847
|
||||
2 0.0 0.083015 0.086807 0.084909 0.084904 0.084702 0.076955 0.076338
|
||||
3 0.0 0.080297 0.084909 0.084414 0.084862 0.085018 0.077793 0.077402
|
||||
4 0.0 0.079157 0.084904 0.084862 0.086076 0.086872 0.079281 0.079383
|
||||
5 0.0 0.078040 0.084702 0.085018 0.086872 0.088212 0.080319 0.080847
|
||||
6 0.0 0.071112 0.076955 0.077793 0.079281 0.080319 0.073719 0.074034
|
||||
7 0.0 0.069847 0.076338 0.077402 0.079383 0.080847 0.074034 0.074693
|
||||
8 0.0 0.068735 0.075760 0.077003 0.079405 0.081233 0.074237 0.075195
|
||||
9 0.0 0.067769 0.075241 0.076628 0.079389 0.081533 0.074375 0.075594
|
||||
10 0.0 0.062136 0.068318 0.069888 0.071921 0.073449 0.067622 0.068376
|
||||
11 0.0 0.061149 0.067723 0.069408 0.071767 0.073583 0.067612 0.068608
|
||||
12 0.0 0.060309 0.067213 0.068991 0.071632 0.073699 0.067599 0.068806
|
||||
13 0.0 0.059601 0.066789 0.068642 0.071529 0.073816 0.067597 0.068992
|
||||
14 0.0 0.059013 0.066449 0.068363 0.071467 0.073949 0.067620 0.069181
|
||||
1 0.0 0.084846 0.071547 0.086679 0.078725 0.071480 0.079609 0.073320
|
||||
2 0.0 0.071547 0.061716 0.073647 0.067908 0.062640 0.068417 0.063857
|
||||
3 0.0 0.086679 0.073647 0.094619 0.086262 0.078697 0.090356 0.083579
|
||||
4 0.0 0.078725 0.067908 0.086262 0.079483 0.073303 0.082874 0.077375
|
||||
5 0.0 0.071480 0.062640 0.078697 0.073303 0.068345 0.076121 0.071741
|
||||
6 0.0 0.079609 0.068417 0.090356 0.082874 0.076121 0.088460 0.082267
|
||||
7 0.0 0.073320 0.063857 0.083579 0.077375 0.071741 0.082267 0.077131
|
||||
8 0.0 0.067787 0.059824 0.077620 0.072521 0.067856 0.076821 0.072597
|
||||
9 0.0 0.062888 0.056235 0.072352 0.068210 0.064388 0.072008 0.068573
|
||||
10 0.0 0.071906 0.062551 0.083694 0.077291 0.071515 0.083361 0.077978
|
||||
11 0.0 0.066654 0.058715 0.077948 0.072611 0.067767 0.078042 0.073548
|
||||
12 0.0 0.062061 0.055344 0.072918 0.068500 0.064461 0.073380 0.069654
|
||||
13 0.0 0.058033 0.052372 0.068505 0.064880 0.061538 0.069285 0.066223
|
||||
14 0.0 0.054491 0.049747 0.064623 0.061685 0.058948 0.065680 0.063192
|
||||
|
||||
8 9 10 11 12 13 14
|
||||
0 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
|
||||
1 0.068735 0.067769 0.062136 0.061149 0.060309 0.059601 0.059013
|
||||
2 0.075760 0.075241 0.068318 0.067723 0.067213 0.066789 0.066449
|
||||
3 0.077003 0.076628 0.069888 0.069408 0.068991 0.068642 0.068363
|
||||
4 0.079405 0.079389 0.071921 0.071767 0.071632 0.071529 0.071467
|
||||
5 0.081233 0.081533 0.073449 0.073583 0.073699 0.073816 0.073949
|
||||
6 0.074237 0.074375 0.067622 0.067612 0.067599 0.067597 0.067620
|
||||
7 0.075195 0.075594 0.068376 0.068608 0.068806 0.068992 0.069181
|
||||
8 0.075959 0.076589 0.068965 0.069409 0.069796 0.070150 0.070491
|
||||
9 0.076589 0.077425 0.069441 0.070074 0.070630 0.071136 0.071614
|
||||
10 0.068965 0.069441 0.063052 0.063364 0.063631 0.063874 0.064110
|
||||
11 0.069409 0.070074 0.063364 0.063851 0.064274 0.064658 0.065020
|
||||
12 0.069796 0.070630 0.063631 0.064274 0.064838 0.065348 0.065826
|
||||
13 0.070150 0.071136 0.063874 0.064658 0.065348 0.065974 0.066558
|
||||
14 0.070491 0.071614 0.064110 0.065020 0.065826 0.066558 0.067240
|
||||
1 0.067787 0.062888 0.071906 0.066654 0.062061 0.058033 0.054491
|
||||
2 0.059824 0.056235 0.062551 0.058715 0.055344 0.052372 0.049747
|
||||
3 0.077620 0.072352 0.083694 0.077948 0.072918 0.068505 0.064623
|
||||
4 0.072521 0.068210 0.077291 0.072611 0.068500 0.064880 0.061685
|
||||
5 0.067856 0.064388 0.071515 0.067767 0.064461 0.061538 0.058948
|
||||
6 0.076821 0.072008 0.083361 0.078042 0.073380 0.069285 0.065680
|
||||
7 0.072597 0.068573 0.077978 0.073548 0.069654 0.066223 0.063192
|
||||
8 0.068852 0.065514 0.073238 0.069578 0.066348 0.063493 0.060963
|
||||
9 0.065514 0.062773 0.069044 0.066051 0.063401 0.061048 0.058955
|
||||
10 0.073238 0.069044 0.079558 0.074881 0.070775 0.067162 0.063977
|
||||
11 0.069578 0.066051 0.074881 0.070959 0.067506 0.064459 0.061765
|
||||
12 0.066348 0.063401 0.070775 0.067506 0.064618 0.062062 0.059795
|
||||
13 0.063493 0.061048 0.067162 0.064459 0.062062 0.059933 0.058038
|
||||
14 0.060963 0.058955 0.063977 0.061765 0.059795 0.058038 0.056468
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
@@ -1439,7 +1435,9 @@ x_{01}x_{00}+x_{11}x_{10} & x_{01}^2+x_{11}^2\\
|
||||
\mathrm{cov}[\boldsymbol{x}_1,\boldsymbol{x}_0] & \mathrm{var}[\boldsymbol{x}_1] \\
|
||||
\end{bmatrix},
|
||||
\end{split}\]</div>
|
||||
<p>where we wrote $<span class="math notranslate nohighlight">\(\boldsymbol{C}[\boldsymbol{x}_0,\boldsymbol{x}_1] = \boldsymbol{C}[\boldsymbol{x}]\)</span><span class="math notranslate nohighlight">\( to indicate that this is the covariance of the vectors \)</span>\boldsymbol{x}<span class="math notranslate nohighlight">\( of the design/feature matrix \)</span>\boldsymbol{X}$.</p>
|
||||
<p>where we wrote <span class="math notranslate nohighlight">\(\boldsymbol{C}[\boldsymbol{x}_0,\boldsymbol{x}_1]=\boldsymbol{C}[\boldsymbol{x}]\)</span> to indicate
|
||||
that this is the covariance of the vectors <span class="math notranslate nohighlight">\(\boldsymbol{x}\)</span> of the
|
||||
design/feature matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\)</span>.</p>
|
||||
<p>It is easy to generalize this to a matrix <span class="math notranslate nohighlight">\(\boldsymbol{X}\in {\mathbb{R}}^{n\times p}\)</span>.</p>
|
||||
</div>
|
||||
<div class="section" id="linking-with-the-svd">
|
||||
@@ -1880,7 +1878,7 @@ Training MSE for OLS
|
||||
3.0
|
||||
</pre></div>
|
||||
</div>
|
||||
<img alt="_images/chapter2_245_1.png" src="_images/chapter2_245_1.png" />
|
||||
<img alt="_images/chapter2_254_1.png" src="_images/chapter2_254_1.png" />
|
||||
</div>
|
||||
</div>
|
||||
<p>We see here that we reach a plateau for the Ridge results. Writing out the coefficients <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span>, we that they are getting smaller and smaller and our error stabilizes since the predicted values of <span class="math notranslate nohighlight">\(\tilde{\boldsymbol{y}}\)</span> approach zero.</p>
|
||||
@@ -2156,7 +2154,7 @@ Training MSE for OLS
|
||||
[ 0. -0.]
|
||||
</pre></div>
|
||||
</div>
|
||||
<img alt="_images/chapter2_247_1.png" src="_images/chapter2_247_1.png" />
|
||||
<img alt="_images/chapter2_256_1.png" src="_images/chapter2_256_1.png" />
|
||||
</div>
|
||||
</div>
|
||||
<p>We bring then back our exponential function example and study all
|
||||
@@ -2259,305 +2257,14 @@ Test MSE OLS
|
||||
0.008675369724975977
|
||||
</pre></div>
|
||||
</div>
|
||||
<img alt="_images/chapter2_249_1.png" src="_images/chapter2_249_1.png" />
|
||||
<img alt="_images/chapter2_258_1.png" src="_images/chapter2_258_1.png" />
|
||||
</div>
|
||||
</div>
|
||||
<p>Both these example send a clear message. The addition of a
|
||||
shrinkage/regularization term implies that we need to perform a search
|
||||
for the optimal values of <span class="math notranslate nohighlight">\(\lambda\)</span>. We will see this throughout these
|
||||
series of lectures.</p>
|
||||
<p>As a small addendum, we note that you can also solve this problem using the convex optimization package <a class="reference external" href="https://cvxopt.org/examples/mlbook/l1regls.html">CVXOPT</a>. This requires, in addition to having installed <strong>CVXOPT</strong>, you need to download the file <em><a class="reference external" href="http://l1regl.py">l1regl.py</a></em>.
|
||||
The following code example solves the simpler problem we discussed above, where we have added the latter python file.</p>
|
||||
<div class="cell docutils container">
|
||||
<div class="cell_input docutils container">
|
||||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">spdiag</span><span class="p">,</span> <span class="n">mul</span><span class="p">,</span> <span class="n">div</span><span class="p">,</span> <span class="n">sqrt</span><span class="p">,</span> <span class="n">normal</span><span class="p">,</span> <span class="n">setseed</span>
|
||||
<span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">blas</span><span class="p">,</span> <span class="n">lapack</span><span class="p">,</span> <span class="n">solvers</span><span class="p">,</span> <span class="n">sparse</span><span class="p">,</span> <span class="n">spmatrix</span>
|
||||
<span class="kn">import</span> <span class="nn">math</span>
|
||||
|
||||
<span class="k">try</span><span class="p">:</span>
|
||||
<span class="kn">import</span> <span class="nn">mosek</span>
|
||||
<span class="kn">import</span> <span class="nn">sys</span>
|
||||
<span class="n">__MOSEK</span> <span class="o">=</span> <span class="kc">True</span>
|
||||
<span class="k">except</span><span class="p">:</span> <span class="n">__MOSEK</span> <span class="o">=</span> <span class="kc">False</span>
|
||||
|
||||
<span class="k">if</span> <span class="n">__MOSEK</span><span class="p">:</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">l1regls_mosek</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
|
||||
<span class="sd">"""</span>
|
||||
|
||||
<span class="sd"> Returns the solution of l1-norm regularized least-squares problem</span>
|
||||
|
||||
<span class="sd"> minimize || A*x - b ||_2^2 + e'*u</span>
|
||||
|
||||
<span class="sd"> subject to -u <= x <= u</span>
|
||||
|
||||
<span class="sd"> """</span>
|
||||
|
||||
<span class="n">m</span><span class="p">,</span> <span class="n">n</span> <span class="o">=</span> <span class="n">A</span><span class="o">.</span><span class="n">size</span>
|
||||
|
||||
<span class="n">env</span> <span class="o">=</span> <span class="n">mosek</span><span class="o">.</span><span class="n">Env</span><span class="p">()</span>
|
||||
<span class="n">task</span> <span class="o">=</span> <span class="n">env</span><span class="o">.</span><span class="n">Task</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">)</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">set_Stream</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">streamtype</span><span class="o">.</span><span class="n">log</span><span class="p">,</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">sys</span><span class="o">.</span><span class="n">stdout</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
|
||||
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">appendvars</span><span class="p">(</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">)</span> <span class="c1"># number of variables</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">appendcons</span><span class="p">(</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">)</span> <span class="c1"># number of constraints</span>
|
||||
|
||||
<span class="c1"># input quadratic objective</span>
|
||||
<span class="n">Q</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="p">(</span><span class="n">n</span><span class="p">,</span><span class="n">n</span><span class="p">))</span>
|
||||
<span class="n">blas</span><span class="o">.</span><span class="n">syrk</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">Q</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">2.0</span><span class="p">,</span> <span class="n">trans</span><span class="o">=</span><span class="s1">'T'</span><span class="p">)</span>
|
||||
|
||||
<span class="n">I</span> <span class="o">=</span> <span class="p">[]</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
|
||||
<span class="n">I</span><span class="o">.</span><span class="n">extend</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="n">i</span><span class="p">,</span><span class="n">n</span><span class="p">))</span>
|
||||
|
||||
<span class="n">J</span> <span class="o">=</span> <span class="p">[]</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
|
||||
<span class="n">J</span><span class="o">.</span><span class="n">extend</span><span class="p">((</span><span class="n">n</span><span class="o">-</span><span class="n">i</span><span class="p">)</span><span class="o">*</span><span class="p">[</span><span class="n">i</span><span class="p">])</span>
|
||||
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putqobj</span><span class="p">(</span><span class="n">I</span><span class="p">,</span> <span class="n">J</span><span class="p">,</span> <span class="nb">list</span><span class="p">(</span><span class="n">Q</span><span class="p">[</span><span class="n">matrix</span><span class="p">(</span><span class="n">I</span><span class="p">)</span> <span class="o">+</span> <span class="n">matrix</span><span class="p">(</span><span class="n">J</span><span class="p">)</span><span class="o">*</span><span class="n">n</span><span class="p">]))</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putclist</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">),</span> <span class="nb">list</span><span class="p">(</span><span class="o">-</span><span class="mi">2</span><span class="o">*</span><span class="n">A</span><span class="o">.</span><span class="n">T</span><span class="o">*</span><span class="n">b</span><span class="p">)</span> <span class="o">+</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">1.0</span><span class="p">])</span> <span class="c1"># setup linear objective</span>
|
||||
|
||||
<span class="c1"># input constraint matrix row by row</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putarow</span><span class="p">(</span> <span class="n">i</span><span class="p">,</span> <span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">],</span> <span class="p">[</span><span class="mf">1.0</span><span class="p">,</span> <span class="o">-</span><span class="mf">1.0</span><span class="p">])</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putarow</span><span class="p">(</span> <span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">,</span> <span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">],</span> <span class="p">[</span><span class="mf">1.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">])</span>
|
||||
|
||||
<span class="c1"># setup bounds on constraints</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putboundslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">accmode</span><span class="o">.</span><span class="n">con</span><span class="p">,</span>
|
||||
<span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="n">mosek</span><span class="o">.</span><span class="n">boundkey</span><span class="o">.</span><span class="n">up</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">])</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putboundslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">accmode</span><span class="o">.</span><span class="n">con</span><span class="p">,</span>
|
||||
<span class="n">n</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="n">mosek</span><span class="o">.</span><span class="n">boundkey</span><span class="o">.</span><span class="n">lo</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">])</span>
|
||||
|
||||
<span class="c1"># setup variable bounds</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putboundslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">accmode</span><span class="o">.</span><span class="n">var</span><span class="p">,</span>
|
||||
<span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="n">mosek</span><span class="o">.</span><span class="n">boundkey</span><span class="o">.</span><span class="n">fr</span><span class="p">],</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">],</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">])</span>
|
||||
|
||||
<span class="c1"># optimize the task</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putobjsense</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">objsense</span><span class="o">.</span><span class="n">minimize</span><span class="p">)</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">optimize</span><span class="p">()</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">solutionsummary</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">streamtype</span><span class="o">.</span><span class="n">log</span><span class="p">)</span>
|
||||
<span class="n">x</span> <span class="o">=</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">]</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">getsolutionslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">soltype</span><span class="o">.</span><span class="n">itr</span><span class="p">,</span> <span class="n">mosek</span><span class="o">.</span><span class="n">solitem</span><span class="o">.</span><span class="n">xx</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
|
||||
|
||||
<span class="k">return</span> <span class="n">matrix</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">l1regls_mosek2</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
|
||||
<span class="sd">"""</span>
|
||||
|
||||
<span class="sd"> Returns the solution of l1-norm regularized least-squares problem</span>
|
||||
|
||||
<span class="sd"> minimize w'*w + e'*u</span>
|
||||
|
||||
<span class="sd"> subject to -u <= x <= u</span>
|
||||
|
||||
<span class="sd"> A*x - w = b</span>
|
||||
|
||||
<span class="sd"> """</span>
|
||||
|
||||
<span class="n">m</span><span class="p">,</span> <span class="n">n</span> <span class="o">=</span> <span class="n">A</span><span class="o">.</span><span class="n">size</span>
|
||||
|
||||
<span class="n">env</span> <span class="o">=</span> <span class="n">mosek</span><span class="o">.</span><span class="n">Env</span><span class="p">()</span>
|
||||
<span class="n">task</span> <span class="o">=</span> <span class="n">env</span><span class="o">.</span><span class="n">Task</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span><span class="mi">0</span><span class="p">)</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">set_Stream</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">streamtype</span><span class="o">.</span><span class="n">log</span><span class="p">,</span> <span class="k">lambda</span> <span class="n">x</span><span class="p">:</span> <span class="n">sys</span><span class="o">.</span><span class="n">stdout</span><span class="o">.</span><span class="n">write</span><span class="p">(</span><span class="n">x</span><span class="p">))</span>
|
||||
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">appendvars</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span> <span class="o">+</span> <span class="n">m</span><span class="p">)</span> <span class="c1"># number of variables</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">appendcons</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span> <span class="o">+</span> <span class="n">m</span><span class="p">)</span> <span class="c1"># number of constraints</span>
|
||||
|
||||
<span class="c1"># input quadratic objective</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putqobj</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">),</span> <span class="nb">range</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">),</span> <span class="n">m</span><span class="o">*</span><span class="p">[</span><span class="mf">2.0</span><span class="p">])</span>
|
||||
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putclist</span><span class="p">(</span><span class="nb">range</span><span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">),</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">]</span> <span class="o">+</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">1.0</span><span class="p">]</span> <span class="o">+</span> <span class="n">m</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">])</span> <span class="c1"># setup linear objective</span>
|
||||
|
||||
<span class="c1"># input constraint matrix row by row</span>
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">):</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putarow</span><span class="p">(</span> <span class="n">i</span><span class="p">,</span> <span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">],</span> <span class="p">[</span><span class="mf">1.0</span><span class="p">,</span> <span class="o">-</span><span class="mf">1.0</span><span class="p">])</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putarow</span><span class="p">(</span> <span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">,</span> <span class="p">[</span><span class="n">i</span><span class="p">,</span> <span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">],</span> <span class="p">[</span><span class="mf">1.0</span><span class="p">,</span> <span class="mf">1.0</span><span class="p">])</span>
|
||||
|
||||
<span class="k">for</span> <span class="n">i</span> <span class="ow">in</span> <span class="nb">range</span><span class="p">(</span><span class="n">m</span><span class="p">):</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putarow</span><span class="p">(</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">,</span> <span class="nb">range</span><span class="p">(</span><span class="n">n</span><span class="p">)</span> <span class="o">+</span> <span class="p">[</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">i</span><span class="p">],</span> <span class="nb">list</span><span class="p">(</span><span class="n">A</span><span class="p">[</span><span class="n">i</span><span class="p">,:])</span> <span class="o">+</span> <span class="p">[</span><span class="o">-</span><span class="mf">1.0</span><span class="p">])</span>
|
||||
|
||||
<span class="c1"># setup bounds on constraints</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putboundslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">accmode</span><span class="o">.</span><span class="n">con</span><span class="p">,</span>
|
||||
<span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="n">mosek</span><span class="o">.</span><span class="n">boundkey</span><span class="o">.</span><span class="n">up</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">])</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putboundslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">accmode</span><span class="o">.</span><span class="n">con</span><span class="p">,</span>
|
||||
<span class="n">n</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="n">mosek</span><span class="o">.</span><span class="n">boundkey</span><span class="o">.</span><span class="n">lo</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">],</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">])</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putboundslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">accmode</span><span class="o">.</span><span class="n">con</span><span class="p">,</span>
|
||||
<span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">,</span> <span class="n">m</span><span class="o">*</span><span class="p">[</span><span class="n">mosek</span><span class="o">.</span><span class="n">boundkey</span><span class="o">.</span><span class="n">fx</span><span class="p">],</span> <span class="nb">list</span><span class="p">(</span><span class="n">b</span><span class="p">),</span> <span class="nb">list</span><span class="p">(</span><span class="n">b</span><span class="p">))</span>
|
||||
|
||||
<span class="c1"># setup variable bounds</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putboundslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">accmode</span><span class="o">.</span><span class="n">var</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">,</span> <span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">)</span><span class="o">*</span><span class="p">[</span><span class="n">mosek</span><span class="o">.</span><span class="n">boundkey</span><span class="o">.</span><span class="n">fr</span><span class="p">],</span>
|
||||
<span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">)</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">],</span> <span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="o">+</span><span class="n">m</span><span class="p">)</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">])</span>
|
||||
|
||||
<span class="c1"># optimize the task</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">putobjsense</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">objsense</span><span class="o">.</span><span class="n">minimize</span><span class="p">)</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">optimize</span><span class="p">()</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">solutionsummary</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">streamtype</span><span class="o">.</span><span class="n">log</span><span class="p">)</span>
|
||||
<span class="n">x</span> <span class="o">=</span> <span class="n">n</span><span class="o">*</span><span class="p">[</span><span class="mf">0.0</span><span class="p">]</span>
|
||||
<span class="n">task</span><span class="o">.</span><span class="n">getsolutionslice</span><span class="p">(</span><span class="n">mosek</span><span class="o">.</span><span class="n">soltype</span><span class="o">.</span><span class="n">itr</span><span class="p">,</span> <span class="n">mosek</span><span class="o">.</span><span class="n">solitem</span><span class="o">.</span><span class="n">xx</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="n">n</span><span class="p">,</span> <span class="n">x</span><span class="p">)</span>
|
||||
|
||||
<span class="k">return</span> <span class="n">matrix</span><span class="p">(</span><span class="n">x</span><span class="p">)</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">l1regls</span><span class="p">(</span><span class="n">A</span><span class="p">,</span> <span class="n">b</span><span class="p">):</span>
|
||||
<span class="sd">"""</span>
|
||||
<span class="sd"> </span>
|
||||
<span class="sd"> Returns the solution of l1-norm regularized least-squares problem</span>
|
||||
<span class="sd"> </span>
|
||||
<span class="sd"> minimize || A*x - b ||_2^2 + || x ||_1.</span>
|
||||
|
||||
<span class="sd"> """</span>
|
||||
|
||||
<span class="n">m</span><span class="p">,</span> <span class="n">n</span> <span class="o">=</span> <span class="n">A</span><span class="o">.</span><span class="n">size</span>
|
||||
<span class="n">q</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">1.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
|
||||
<span class="n">q</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="o">-</span><span class="mf">2.0</span> <span class="o">*</span> <span class="n">A</span><span class="o">.</span><span class="n">T</span> <span class="o">*</span> <span class="n">b</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">P</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">v</span><span class="p">,</span> <span class="n">alpha</span> <span class="o">=</span> <span class="mf">1.0</span><span class="p">,</span> <span class="n">beta</span> <span class="o">=</span> <span class="mf">0.0</span> <span class="p">):</span>
|
||||
<span class="sd">"""</span>
|
||||
<span class="sd"> v := alpha * 2.0 * [ A'*A, 0; 0, 0 ] * u + beta * v </span>
|
||||
<span class="sd"> """</span>
|
||||
<span class="n">v</span> <span class="o">*=</span> <span class="n">beta</span>
|
||||
<span class="n">v</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">+=</span> <span class="n">alpha</span> <span class="o">*</span> <span class="mf">2.0</span> <span class="o">*</span> <span class="n">A</span><span class="o">.</span><span class="n">T</span> <span class="o">*</span> <span class="p">(</span><span class="n">A</span> <span class="o">*</span> <span class="n">u</span><span class="p">[:</span><span class="n">n</span><span class="p">])</span>
|
||||
|
||||
|
||||
<span class="k">def</span> <span class="nf">G</span><span class="p">(</span><span class="n">u</span><span class="p">,</span> <span class="n">v</span><span class="p">,</span> <span class="n">alpha</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">beta</span><span class="o">=</span><span class="mf">0.0</span><span class="p">,</span> <span class="n">trans</span><span class="o">=</span><span class="s1">'N'</span><span class="p">):</span>
|
||||
<span class="sd">"""</span>
|
||||
<span class="sd"> v := alpha*[I, -I; -I, -I] * u + beta * v (trans = 'N' or 'T')</span>
|
||||
<span class="sd"> """</span>
|
||||
|
||||
<span class="n">v</span> <span class="o">*=</span> <span class="n">beta</span>
|
||||
<span class="n">v</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">+=</span> <span class="n">alpha</span><span class="o">*</span><span class="p">(</span><span class="n">u</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">-</span> <span class="n">u</span><span class="p">[</span><span class="n">n</span><span class="p">:])</span>
|
||||
<span class="n">v</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="o">+=</span> <span class="n">alpha</span><span class="o">*</span><span class="p">(</span><span class="o">-</span><span class="n">u</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">-</span> <span class="n">u</span><span class="p">[</span><span class="n">n</span><span class="p">:])</span>
|
||||
|
||||
<span class="n">h</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="p">(</span><span class="mi">2</span><span class="o">*</span><span class="n">n</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
|
||||
|
||||
|
||||
<span class="c1"># Customized solver for the KKT system </span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># [ 2.0*A'*A 0 I -I ] [x[:n] ] [bx[:n] ]</span>
|
||||
<span class="c1"># [ 0 0 -I -I ] [x[n:] ] = [bx[n:] ].</span>
|
||||
<span class="c1"># [ I -I -D1^-1 0 ] [zl[:n]] [bzl[:n]]</span>
|
||||
<span class="c1"># [ -I -I 0 -D2^-1 ] [zl[n:]] [bzl[n:]]</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># where D1 = W['di'][:n]**2, D2 = W['di'][:n]**2.</span>
|
||||
<span class="c1"># </span>
|
||||
<span class="c1"># We first eliminate zl and x[n:]:</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># ( 2*A'*A + 4*D1*D2*(D1+D2)^-1 ) * x[:n] = </span>
|
||||
<span class="c1"># bx[:n] - (D2-D1)*(D1+D2)^-1 * bx[n:] + </span>
|
||||
<span class="c1"># D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] - </span>
|
||||
<span class="c1"># D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:] </span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># x[n:] = (D1+D2)^-1 * ( bx[n:] - D1*bzl[:n] - D2*bzl[n:] ) </span>
|
||||
<span class="c1"># - (D2-D1)*(D1+D2)^-1 * x[:n] </span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># zl[:n] = D1 * ( x[:n] - x[n:] - bzl[:n] )</span>
|
||||
<span class="c1"># zl[n:] = D2 * (-x[:n] - x[n:] - bzl[n:] ).</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># The first equation has the form</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># (A'*A + D)*x[:n] = rhs</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># and is equivalent to</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># [ D A' ] [ x:n] ] = [ rhs ]</span>
|
||||
<span class="c1"># [ A -I ] [ v ] [ 0 ].</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># It can be solved as </span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># ( A*D^-1*A' + I ) * v = A * D^-1 * rhs</span>
|
||||
<span class="c1"># x[:n] = D^-1 * ( rhs - A'*v ).</span>
|
||||
|
||||
<span class="n">S</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="p">(</span><span class="n">m</span><span class="p">,</span><span class="n">m</span><span class="p">))</span>
|
||||
<span class="n">Asc</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="p">(</span><span class="n">m</span><span class="p">,</span><span class="n">n</span><span class="p">))</span>
|
||||
<span class="n">v</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span><span class="mf">0.0</span><span class="p">,</span> <span class="p">(</span><span class="n">m</span><span class="p">,</span><span class="mi">1</span><span class="p">))</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">Fkkt</span><span class="p">(</span><span class="n">W</span><span class="p">):</span>
|
||||
|
||||
<span class="c1"># Factor </span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># S = A*D^-1*A' + I </span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># where D = 2*D1*D2*(D1+D2)^-1, D1 = d[:n]**-2, D2 = d[n:]**-2.</span>
|
||||
|
||||
<span class="n">d1</span><span class="p">,</span> <span class="n">d2</span> <span class="o">=</span> <span class="n">W</span><span class="p">[</span><span class="s1">'di'</span><span class="p">][:</span><span class="n">n</span><span class="p">]</span><span class="o">**</span><span class="mi">2</span><span class="p">,</span> <span class="n">W</span><span class="p">[</span><span class="s1">'di'</span><span class="p">][</span><span class="n">n</span><span class="p">:]</span><span class="o">**</span><span class="mi">2</span>
|
||||
|
||||
<span class="c1"># ds is square root of diagonal of D</span>
|
||||
<span class="n">ds</span> <span class="o">=</span> <span class="n">math</span><span class="o">.</span><span class="n">sqrt</span><span class="p">(</span><span class="mf">2.0</span><span class="p">)</span> <span class="o">*</span> <span class="n">div</span><span class="p">(</span> <span class="n">mul</span><span class="p">(</span> <span class="n">W</span><span class="p">[</span><span class="s1">'di'</span><span class="p">][:</span><span class="n">n</span><span class="p">],</span> <span class="n">W</span><span class="p">[</span><span class="s1">'di'</span><span class="p">][</span><span class="n">n</span><span class="p">:]),</span>
|
||||
<span class="n">sqrt</span><span class="p">(</span><span class="n">d1</span><span class="o">+</span><span class="n">d2</span><span class="p">)</span> <span class="p">)</span>
|
||||
<span class="n">d3</span> <span class="o">=</span> <span class="n">div</span><span class="p">(</span><span class="n">d2</span> <span class="o">-</span> <span class="n">d1</span><span class="p">,</span> <span class="n">d1</span> <span class="o">+</span> <span class="n">d2</span><span class="p">)</span>
|
||||
|
||||
<span class="c1"># Asc = A*diag(d)^-1/2</span>
|
||||
<span class="n">Asc</span> <span class="o">=</span> <span class="n">A</span> <span class="o">*</span> <span class="n">spdiag</span><span class="p">(</span><span class="n">ds</span><span class="o">**-</span><span class="mi">1</span><span class="p">)</span>
|
||||
|
||||
<span class="c1"># S = I + A * D^-1 * A'</span>
|
||||
<span class="n">blas</span><span class="o">.</span><span class="n">syrk</span><span class="p">(</span><span class="n">Asc</span><span class="p">,</span> <span class="n">S</span><span class="p">)</span>
|
||||
<span class="n">S</span><span class="p">[::</span><span class="n">m</span><span class="o">+</span><span class="mi">1</span><span class="p">]</span> <span class="o">+=</span> <span class="mf">1.0</span>
|
||||
<span class="n">lapack</span><span class="o">.</span><span class="n">potrf</span><span class="p">(</span><span class="n">S</span><span class="p">)</span>
|
||||
|
||||
<span class="k">def</span> <span class="nf">g</span><span class="p">(</span><span class="n">x</span><span class="p">,</span> <span class="n">y</span><span class="p">,</span> <span class="n">z</span><span class="p">):</span>
|
||||
|
||||
<span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="mf">0.5</span> <span class="o">*</span> <span class="p">(</span> <span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">-</span> <span class="n">mul</span><span class="p">(</span><span class="n">d3</span><span class="p">,</span> <span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">:])</span> <span class="o">+</span>
|
||||
<span class="n">mul</span><span class="p">(</span><span class="n">d1</span><span class="p">,</span> <span class="n">z</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">+</span> <span class="n">mul</span><span class="p">(</span><span class="n">d3</span><span class="p">,</span> <span class="n">z</span><span class="p">[:</span><span class="n">n</span><span class="p">]))</span> <span class="o">-</span> <span class="n">mul</span><span class="p">(</span><span class="n">d2</span><span class="p">,</span> <span class="n">z</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="o">-</span>
|
||||
<span class="n">mul</span><span class="p">(</span><span class="n">d3</span><span class="p">,</span> <span class="n">z</span><span class="p">[</span><span class="n">n</span><span class="p">:]))</span> <span class="p">)</span>
|
||||
<span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="n">div</span><span class="p">(</span> <span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">],</span> <span class="n">ds</span><span class="p">)</span>
|
||||
|
||||
<span class="c1"># Solve</span>
|
||||
<span class="c1">#</span>
|
||||
<span class="c1"># S * v = 0.5 * A * D^-1 * ( bx[:n] - </span>
|
||||
<span class="c1"># (D2-D1)*(D1+D2)^-1 * bx[n:] + </span>
|
||||
<span class="c1"># D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] - </span>
|
||||
<span class="c1"># D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:] )</span>
|
||||
|
||||
<span class="n">blas</span><span class="o">.</span><span class="n">gemv</span><span class="p">(</span><span class="n">Asc</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">v</span><span class="p">)</span>
|
||||
<span class="n">lapack</span><span class="o">.</span><span class="n">potrs</span><span class="p">(</span><span class="n">S</span><span class="p">,</span> <span class="n">v</span><span class="p">)</span>
|
||||
|
||||
<span class="c1"># x[:n] = D^-1 * ( rhs - A'*v ).</span>
|
||||
<span class="n">blas</span><span class="o">.</span><span class="n">gemv</span><span class="p">(</span><span class="n">Asc</span><span class="p">,</span> <span class="n">v</span><span class="p">,</span> <span class="n">x</span><span class="p">,</span> <span class="n">alpha</span><span class="o">=-</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">beta</span><span class="o">=</span><span class="mf">1.0</span><span class="p">,</span> <span class="n">trans</span><span class="o">=</span><span class="s1">'T'</span><span class="p">)</span>
|
||||
<span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="n">div</span><span class="p">(</span><span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">],</span> <span class="n">ds</span><span class="p">)</span>
|
||||
|
||||
<span class="c1"># x[n:] = (D1+D2)^-1 * ( bx[n:] - D1*bzl[:n] - D2*bzl[n:] ) </span>
|
||||
<span class="c1"># - (D2-D1)*(D1+D2)^-1 * x[:n] </span>
|
||||
<span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="o">=</span> <span class="n">div</span><span class="p">(</span> <span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="o">-</span> <span class="n">mul</span><span class="p">(</span><span class="n">d1</span><span class="p">,</span> <span class="n">z</span><span class="p">[:</span><span class="n">n</span><span class="p">])</span> <span class="o">-</span> <span class="n">mul</span><span class="p">(</span><span class="n">d2</span><span class="p">,</span> <span class="n">z</span><span class="p">[</span><span class="n">n</span><span class="p">:]),</span> <span class="n">d1</span><span class="o">+</span><span class="n">d2</span> <span class="p">)</span>\
|
||||
<span class="o">-</span> <span class="n">mul</span><span class="p">(</span> <span class="n">d3</span><span class="p">,</span> <span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="p">)</span>
|
||||
|
||||
<span class="c1"># zl[:n] = D1^1/2 * ( x[:n] - x[n:] - bzl[:n] )</span>
|
||||
<span class="c1"># zl[n:] = D2^1/2 * ( -x[:n] - x[n:] - bzl[n:] ).</span>
|
||||
<span class="n">z</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">=</span> <span class="n">mul</span><span class="p">(</span> <span class="n">W</span><span class="p">[</span><span class="s1">'di'</span><span class="p">][:</span><span class="n">n</span><span class="p">],</span> <span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">-</span> <span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="o">-</span> <span class="n">z</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="p">)</span>
|
||||
<span class="n">z</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="o">=</span> <span class="n">mul</span><span class="p">(</span> <span class="n">W</span><span class="p">[</span><span class="s1">'di'</span><span class="p">][</span><span class="n">n</span><span class="p">:],</span> <span class="o">-</span><span class="n">x</span><span class="p">[:</span><span class="n">n</span><span class="p">]</span> <span class="o">-</span> <span class="n">x</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="o">-</span> <span class="n">z</span><span class="p">[</span><span class="n">n</span><span class="p">:]</span> <span class="p">)</span>
|
||||
|
||||
<span class="k">return</span> <span class="n">g</span>
|
||||
|
||||
<span class="k">return</span> <span class="n">solvers</span><span class="o">.</span><span class="n">coneqp</span><span class="p">(</span><span class="n">P</span><span class="p">,</span> <span class="n">q</span><span class="p">,</span> <span class="n">G</span><span class="p">,</span> <span class="n">h</span><span class="p">,</span> <span class="n">kktsolver</span> <span class="o">=</span> <span class="n">Fkkt</span><span class="p">)[</span><span class="s1">'x'</span><span class="p">][:</span><span class="n">n</span><span class="p">]</span>
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output traceback highlight-ipythontb notranslate"><div class="highlight"><pre><span></span><span class="gt">---------------------------------------------------------------------------</span>
|
||||
<span class="ne">ModuleNotFoundError</span><span class="g g-Whitespace"> </span>Traceback (most recent call last)
|
||||
<span class="o"><</span><span class="n">ipython</span><span class="o">-</span><span class="nb">input</span><span class="o">-</span><span class="mi">12</span><span class="o">-</span><span class="n">d670a873ab0c</span><span class="o">></span> <span class="ow">in</span> <span class="o"><</span><span class="n">module</span><span class="o">></span>
|
||||
<span class="ne">----> </span><span class="mi">1</span> <span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">spdiag</span><span class="p">,</span> <span class="n">mul</span><span class="p">,</span> <span class="n">div</span><span class="p">,</span> <span class="n">sqrt</span><span class="p">,</span> <span class="n">normal</span><span class="p">,</span> <span class="n">setseed</span>
|
||||
<span class="g g-Whitespace"> </span><span class="mi">2</span> <span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">blas</span><span class="p">,</span> <span class="n">lapack</span><span class="p">,</span> <span class="n">solvers</span><span class="p">,</span> <span class="n">sparse</span><span class="p">,</span> <span class="n">spmatrix</span>
|
||||
<span class="g g-Whitespace"> </span><span class="mi">3</span> <span class="kn">import</span> <span class="nn">math</span>
|
||||
<span class="g g-Whitespace"> </span><span class="mi">4</span>
|
||||
<span class="g g-Whitespace"> </span><span class="mi">5</span> <span class="k">try</span><span class="p">:</span>
|
||||
|
||||
<span class="ne">ModuleNotFoundError</span>: No module named 'cvxopt'
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
</div>
|
||||
<p>Then we call the above functions and solve the problem, as done here</p>
|
||||
<div class="cell docutils container">
|
||||
<div class="cell_input docutils container">
|
||||
<div class="highlight-ipython3 notranslate"><div class="highlight"><pre><span></span><span class="kn">from</span> <span class="nn">cvxopt</span> <span class="kn">import</span> <span class="n">matrix</span><span class="p">,</span> <span class="n">normal</span>
|
||||
|
||||
<span class="n">X</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span> <span class="p">[</span> <span class="p">[</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">],</span> <span class="p">[</span><span class="mi">0</span><span class="p">,</span> <span class="mi">1</span><span class="p">,</span> <span class="mi">3</span><span class="p">]])</span>
|
||||
<span class="n">y</span> <span class="o">=</span> <span class="n">matrix</span><span class="p">(</span> <span class="p">[</span><span class="mi">4</span><span class="p">,</span> <span class="mi">2</span><span class="p">,</span> <span class="mi">3</span><span class="p">])</span>
|
||||
<span class="n">x</span> <span class="o">=</span> <span class="n">l1regls</span><span class="p">(</span><span class="n">X</span><span class="p">,</span><span class="n">y</span><span class="p">)</span>
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
</div>
|
||||
<p><strong>More text will be added to this example.</strong></p>
|
||||
<p>As a small addendum, we note that you can also solve this problem using the convex optimization package <a class="reference external" href="https://cvxopt.org/examples/mlbook/l1regls.html">CVXOPT</a>. This requires, in addition to having installed <strong>CVXOPT</strong>, you need to download the file <em><a class="reference external" href="http://l1regl.py">l1regl.py</a></em>.</p>
|
||||
</div>
|
||||
<div class="section" id="linking-the-regression-analysis-with-a-statistical-interpretation">
|
||||
<h2><span class="section-number">4.11. </span>Linking the regression analysis with a statistical interpretation<a class="headerlink" href="#linking-the-regression-analysis-with-a-statistical-interpretation" title="Permalink to this headline">¶</a></h2>
|
||||
@@ -2971,6 +2678,14 @@ order to another one.</p>
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output stream highlight-myst-ansi notranslate"><div class="highlight"><pre><span></span>[ 1.0169643 0.27924636 -1.4087793 1.03308408 0. ]
|
||||
Test MSE OLS
|
||||
0.958228616652075
|
||||
</pre></div>
|
||||
</div>
|
||||
<img alt="_images/chapter2_324_1.png" src="_images/chapter2_324_1.png" />
|
||||
</div>
|
||||
</div>
|
||||
<p>How can we understand this?</p>
|
||||
<p>Let us write out the values of the coefficients <span class="math notranslate nohighlight">\(\beta_i\)</span> as functions
|
||||
@@ -3032,6 +2747,228 @@ large variance (normally for higher orders in the polynomial).</p>
|
||||
</pre></div>
|
||||
</div>
|
||||
</div>
|
||||
<div class="cell_output docutils container">
|
||||
<div class="output text_html"><div>
|
||||
<style scoped>
|
||||
.dataframe tbody tr th:only-of-type {
|
||||
vertical-align: middle;
|
||||
}
|
||||
|
||||
.dataframe tbody tr th {
|
||||
vertical-align: top;
|
||||
}
|
||||
|
||||
.dataframe thead th {
|
||||
text-align: right;
|
||||
}
|
||||
</style>
|
||||
<table border="1" class="dataframe">
|
||||
<thead>
|
||||
<tr style="text-align: right;">
|
||||
<th></th>
|
||||
<th>beta</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<th>0</th>
|
||||
<td>0.986699</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>1</th>
|
||||
<td>-0.606760</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>2</th>
|
||||
<td>1.280573</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>3</th>
|
||||
<td>-0.850164</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>4</th>
|
||||
<td>0.000000</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
</div></div><div class="output text_html"><div>
|
||||
<style scoped>
|
||||
.dataframe tbody tr th:only-of-type {
|
||||
vertical-align: middle;
|
||||
}
|
||||
|
||||
.dataframe tbody tr th {
|
||||
vertical-align: top;
|
||||
}
|
||||
|
||||
.dataframe thead th {
|
||||
text-align: right;
|
||||
}
|
||||
</style>
|
||||
<table border="1" class="dataframe">
|
||||
<thead>
|
||||
<tr style="text-align: right;">
|
||||
<th></th>
|
||||
<th>beta</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<th>0</th>
|
||||
<td>0.978553</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>1</th>
|
||||
<td>-0.511888</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>2</th>
|
||||
<td>1.051418</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>3</th>
|
||||
<td>-0.701370</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>4</th>
|
||||
<td>0.000000</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
</div></div><div class="output text_html"><div>
|
||||
<style scoped>
|
||||
.dataframe tbody tr th:only-of-type {
|
||||
vertical-align: middle;
|
||||
}
|
||||
|
||||
.dataframe tbody tr th {
|
||||
vertical-align: top;
|
||||
}
|
||||
|
||||
.dataframe thead th {
|
||||
text-align: right;
|
||||
}
|
||||
</style>
|
||||
<table border="1" class="dataframe">
|
||||
<thead>
|
||||
<tr style="text-align: right;">
|
||||
<th></th>
|
||||
<th>beta</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<th>0</th>
|
||||
<td>0.946957</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>1</th>
|
||||
<td>-0.162246</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>2</th>
|
||||
<td>0.221921</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>3</th>
|
||||
<td>-0.167787</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>4</th>
|
||||
<td>0.000000</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
</div></div><div class="output text_html"><div>
|
||||
<style scoped>
|
||||
.dataframe tbody tr th:only-of-type {
|
||||
vertical-align: middle;
|
||||
}
|
||||
|
||||
.dataframe tbody tr th {
|
||||
vertical-align: top;
|
||||
}
|
||||
|
||||
.dataframe thead th {
|
||||
text-align: right;
|
||||
}
|
||||
</style>
|
||||
<table border="1" class="dataframe">
|
||||
<thead>
|
||||
<tr style="text-align: right;">
|
||||
<th></th>
|
||||
<th>beta</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<th>0</th>
|
||||
<td>0.906747</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>1</th>
|
||||
<td>0.017665</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>2</th>
|
||||
<td>-0.029483</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>3</th>
|
||||
<td>-0.053849</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>4</th>
|
||||
<td>0.000000</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
</div></div><div class="output text_html"><div>
|
||||
<style scoped>
|
||||
.dataframe tbody tr th:only-of-type {
|
||||
vertical-align: middle;
|
||||
}
|
||||
|
||||
.dataframe tbody tr th {
|
||||
vertical-align: top;
|
||||
}
|
||||
|
||||
.dataframe thead th {
|
||||
text-align: right;
|
||||
}
|
||||
</style>
|
||||
<table border="1" class="dataframe">
|
||||
<thead>
|
||||
<tr style="text-align: right;">
|
||||
<th></th>
|
||||
<th>beta</th>
|
||||
</tr>
|
||||
</thead>
|
||||
<tbody>
|
||||
<tr>
|
||||
<th>0</th>
|
||||
<td>0.718165</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>1</th>
|
||||
<td>0.156956</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>2</th>
|
||||
<td>0.040102</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>3</th>
|
||||
<td>-0.001880</td>
|
||||
</tr>
|
||||
<tr>
|
||||
<th>4</th>
|
||||
<td>0.000000</td>
|
||||
</tr>
|
||||
</tbody>
|
||||
</table>
|
||||
</div></div></div>
|
||||
</div>
|
||||
<p>As an exercise, repeat these calculations with ordinary least squares
|
||||
only with and without noise. Calculate thereafter the variance of the
|
||||
@@ -3074,8 +3011,7 @@ already modeled and an unknown prior, we are now ready to make
|
||||
additional models for the prior.</p>
|
||||
<p>We can, based on our discussions of the variance of <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> and
|
||||
the mean value, assume that the prior for the values <span class="math notranslate nohighlight">\(\boldsymbol{\beta}\)</span> is
|
||||
given by a Gaussian with mean value zero and variance <span class="math notranslate nohighlight">\(\tau^2\)</span>, that
|
||||
is</p>
|
||||
given by a Gaussian with mean value zero and variance <span class="math notranslate nohighlight">\(\tau^2\)</span>, that</p>
|
||||
<div class="math notranslate nohighlight">
|
||||
\[
|
||||
p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}.
|
||||
|
||||
File diff suppressed because one or more lines are too long
File diff suppressed because one or more lines are too long
@@ -1,15 +1,15 @@
|
||||
#!/usr/bin/env python
|
||||
# coding: utf-8
|
||||
|
||||
# <!-- HTML file automatically generated from DocOnce source (https://github.com/doconce/doconce/)
|
||||
# doconce format html chapter2.do.txt -->
|
||||
|
||||
# # Ridge and Lasso Regression
|
||||
#
|
||||
#
|
||||
#
|
||||
|
||||
# ## Mathematical Interpretation of Ordinary Least Squares
|
||||
#
|
||||
# What is presented here is a mathematical analysis of various regression algorithms (ordinary least squares, Ridge and Lasso Regression). The analysis is based on an important algorithm in linear algebra, the so-called Singular Value Decomposition (SVD).
|
||||
#
|
||||
#
|
||||
# We have shown that in ordinary least squares (OLS) the optimal parameters $\beta$ are given by
|
||||
|
||||
# $$
|
||||
@@ -39,9 +39,6 @@
|
||||
# The matrix $\boldsymbol{A}$ has the important property that $\boldsymbol{A}^2=\boldsymbol{A}$. This is the definition of a [projection matrix](https://en.wikipedia.org/wiki/Projection_matrix).
|
||||
# We can then interpret our optimal model $\tilde{\boldsymbol{y}}$ as being represented by an orthogonal projection of $\boldsymbol{y}$ onto a space defined by the column vectors of $\boldsymbol{X}$. In our case here the matrix $\boldsymbol{A}$ is a square matrix. If it is a general rectangular matrix we have an oblique projection matrix.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# We have defined the residual error as
|
||||
|
||||
# $$
|
||||
@@ -50,7 +47,6 @@
|
||||
|
||||
# The residual errors are then the projections of $\boldsymbol{y}$ onto the orthogonal component of the space defined by the column vectors of $\boldsymbol{X}$.
|
||||
#
|
||||
#
|
||||
# If the matrix $\boldsymbol{X}$ is an orthogonal (or unitary in case of complex values) matrix, we have
|
||||
|
||||
# $$
|
||||
@@ -69,14 +65,10 @@
|
||||
# \boldsymbol{\epsilon}=\boldsymbol{y}-\tilde{\boldsymbol{y}}=0.
|
||||
# $$
|
||||
|
||||
# This serves also as a useful test of our codes.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# This serves also as a useful test of our codes.
|
||||
|
||||
# ## The singular value decomposition
|
||||
#
|
||||
#
|
||||
# The examples we have looked at so far are cases where we normally can
|
||||
# invert the matrix $\boldsymbol{X}^T\boldsymbol{X}$. Using a polynomial expansion where we fit of various functions leads to
|
||||
# row vectors of the design matrix which are essentially orthogonal due
|
||||
@@ -84,7 +76,6 @@
|
||||
# design matrix is then often done via a so-called LU, QR or Cholesky
|
||||
# decomposition.
|
||||
#
|
||||
#
|
||||
# As we will also see in the first project,
|
||||
# this may
|
||||
# however not the be case in general and a standard matrix inversion
|
||||
@@ -108,9 +99,6 @@
|
||||
# in the principal component analysis where high-dimensional data can be
|
||||
# reduced to the statistically relevant features.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# One of the typical problems we encounter with linear regression, in particular
|
||||
# when the matrix $\boldsymbol{X}$ (our so-called design matrix) is high-dimensional,
|
||||
# are problems with near singular or singular matrices. The column vectors of $\boldsymbol{X}$
|
||||
@@ -156,8 +144,6 @@
|
||||
# We see easily that $\mbox{det}(\boldsymbol{X}) = x_{11} x_{22} - x_{12} x_{21} = 1 \times (-1) - 1 \times (-1) = 0$. Hence, $\mathbf{X}$ is singular and its inverse is undefined.
|
||||
# This is equivalent to saying that the matrix $\boldsymbol{X}$ has at least an eigenvalue which is zero.
|
||||
#
|
||||
#
|
||||
#
|
||||
# If our design matrix $\boldsymbol{X}$ which enters the linear regression problem
|
||||
|
||||
# <!-- Equation labels as ordinary links -->
|
||||
@@ -182,14 +168,10 @@
|
||||
# \boldsymbol{X}^{T} \boldsymbol{X} \rightarrow \boldsymbol{X}^{T} \boldsymbol{X}+\lambda \boldsymbol{I},
|
||||
# $$
|
||||
|
||||
# where $\boldsymbol{I}$ is the identity matrix. When we discuss **Ridge** regression this is actually what we end up evaluating. The parameter $\lambda$ is called a hyperparameter. More about this later.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# where $\boldsymbol{I}$ is the identity matrix. When we discuss **Ridge** regression this is actually what we end up evaluating. The parameter $\lambda$ is called a hyperparameter. More about this later.
|
||||
|
||||
# ## Basic math of the SVD
|
||||
#
|
||||
#
|
||||
# From standard linear algebra we know that a square matrix $\boldsymbol{X}$ can be diagonalized if and only it is
|
||||
# a so-called [normal matrix](https://en.wikipedia.org/wiki/Normal_matrix), that is if $\boldsymbol{X}\in {\mathbb{R}}^{n\times n}$
|
||||
# we have $\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{X}^T\boldsymbol{X}$ or if $\boldsymbol{X}\in {\mathbb{C}}^{n\times n}$ we have $\boldsymbol{X}\boldsymbol{X}^{\dagger}=\boldsymbol{X}^{\dagger}\boldsymbol{X}$.
|
||||
@@ -225,10 +207,6 @@
|
||||
# is not diagonalizable, it is a so-called [defective matrix](https://en.wikipedia.org/wiki/Defective_matrix). It is easy to see that the condition
|
||||
# $\boldsymbol{X}\boldsymbol{X}^T=\boldsymbol{X}^T\boldsymbol{X}$ is not fulfilled.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# However, and this is the strength of the SVD algorithm, any general
|
||||
# matrix $\boldsymbol{X}$ can be decomposed in terms of a diagonal matrix and
|
||||
# two orthogonal/unitary matrices. The [Singular Value Decompostion
|
||||
@@ -270,7 +248,6 @@
|
||||
#
|
||||
# The columns of $\boldsymbol{U}$ are called the left singular vectors while the columns of $\boldsymbol{V}$ are the right singular vectors.
|
||||
#
|
||||
#
|
||||
# If we assume that $n > p$, then our matrix $\boldsymbol{U}$ has dimension $n
|
||||
# \times n$. The last $n-p$ columns of $\boldsymbol{U}$ become however
|
||||
# irrelevant in our calculations since they are multiplied with the
|
||||
@@ -287,7 +264,7 @@
|
||||
# If $p > n$, then only the first $n$ columns of $\boldsymbol{V}$ are computed and $\boldsymbol{\Sigma}$ has dimension $n\times n$.
|
||||
# The $n=p$ case is obvious, we retain the full SVD.
|
||||
# In general the economy-size SVD leads to less FLOPS and still conserving the desired accuracy.
|
||||
#
|
||||
|
||||
# ## Codes for the SVD
|
||||
|
||||
# In[1]:
|
||||
@@ -333,8 +310,6 @@ print(C-X)
|
||||
# inversion algorithm for matrix inversion with $\boldsymbol{X}^T\boldsymbol{X}$ results
|
||||
# in the program terminating due to a singular matrix.
|
||||
#
|
||||
#
|
||||
#
|
||||
# The $U$, $S$, and $V$ matrices returned from the **svd()** function
|
||||
# cannot be multiplied directly.
|
||||
#
|
||||
@@ -347,8 +322,7 @@ print(C-X)
|
||||
# If you wish to include the zero singular values, you will need to
|
||||
# resize the matrices and set up a diagonal matrix as done in the above
|
||||
# example
|
||||
#
|
||||
#
|
||||
|
||||
# ## Code for SVD and Inversion of Matrices
|
||||
#
|
||||
# How do we use the SVD to invert a matrix $\boldsymbol{X}^T\boldsymbol{X}$ which is singular or near singular?
|
||||
@@ -410,7 +384,7 @@ print(np.abs(B-C))
|
||||
# \boldsymbol{A}_{\mathrm{PI}}= \boldsymbol{V}\boldsymbol{D}_{\mathrm{PI}}\boldsymbol{U}^T,
|
||||
# $$
|
||||
|
||||
# where $\boldsymbol{D}_{\mathrm{PI}}$ can be calculated by creating a diagonal matrix from $\boldsymbol{Sigma}$ where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.
|
||||
# where $\boldsymbol{D}_{\mathrm{PI}}$ can be calculated by creating a diagonal matrix from $\boldsymbol{\Sigma}$ where we only keep the singular values (the non-zero values). The following code computes the pseudoinvers of the matrix based on the SVD.
|
||||
|
||||
# In[4]:
|
||||
|
||||
@@ -441,11 +415,7 @@ print(np.abs(C-B))
|
||||
|
||||
|
||||
# As you can see from these examples, our own decomposition based on the SVD agrees the pseudoinverse algorithm provided by **Numpy**.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
|
||||
# ## Mathematics of the SVD and implications
|
||||
#
|
||||
# Let us take a closer look at the mathematics of the SVD and the various implications for machine learning studies.
|
||||
@@ -481,7 +451,6 @@ print(np.abs(C-B))
|
||||
|
||||
# All values beyond $p-1$ are all zero.
|
||||
#
|
||||
#
|
||||
# As an example, consider the following $3\times 2$ example for the matrix $\boldsymbol{\Sigma}$
|
||||
|
||||
# $$
|
||||
@@ -538,8 +507,6 @@ print(np.abs(C-B))
|
||||
# contain only zeros. This will have important consequences for our SVD
|
||||
# decomposition of the design matrix.
|
||||
#
|
||||
#
|
||||
#
|
||||
# The matrix that may cause problems for us is $\boldsymbol{X}^T\boldsymbol{X}$. Using the SVD we can rewrite this matrix as
|
||||
|
||||
# $$
|
||||
@@ -554,12 +521,6 @@ print(np.abs(C-B))
|
||||
|
||||
# We define $\boldsymbol{\Sigma}^T\boldsymbol{\Sigma}=\tilde{\boldsymbol{\Sigma}}^2$ which is a diagonal matrix containing only the singular values squared. It has dimensionality $p \times p$.
|
||||
#
|
||||
# This means, using the orthogonality of $\boldsymbol{V}$, that we get
|
||||
|
||||
# $$
|
||||
# \boldsymbol{X}^T\boldsymbol{X}=\tilde{\boldsymbol{\Sigma}}^2.
|
||||
# $$
|
||||
|
||||
# We can now insert the result for the matrix $\boldsymbol{X}^T\boldsymbol{X}$ into our equation for ordinary least squares where
|
||||
|
||||
# $$
|
||||
@@ -569,10 +530,10 @@ print(np.abs(C-B))
|
||||
# and using our SVD decomposition of $\boldsymbol{X}$ we have
|
||||
|
||||
# $$
|
||||
# \tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T\tilde{\boldsymbol{\Sigma}}^{-2}\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{y},
|
||||
# \tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{\Sigma}\boldsymbol{V}^T\left(\boldsymbol{V}\tilde{\boldsymbol{\Sigma}}^{2}(\boldsymbol{V}^T\right)^{-1}\boldsymbol{V}\boldsymbol{\Sigma}^T\boldsymbol{U}^T\boldsymbol{y},
|
||||
# $$
|
||||
|
||||
# which gives us, using the orthogonality of the matrices $\boldsymbol{U}$ and $\boldsymbol{V}$,
|
||||
# which gives us, using the orthogonality of the matrices $\boldsymbol{U}$ and $\boldsymbol{V}$,,
|
||||
|
||||
# $$
|
||||
# \tilde{y}_{\mathrm{OLS}}=\boldsymbol{U}\boldsymbol{U}^T\boldsymbol{y}=\sum_{i=0}^{p-1}\boldsymbol{u}_i\boldsymbol{u}^T_j\boldsymbol{y},
|
||||
@@ -587,8 +548,7 @@ print(np.abs(C-B))
|
||||
# that belong to $i>p-1$, result in only zeros when we perform the multiplications. This means that the sum above has non-zero elements only up to $i=p-1$. This corresponds also to the number of singular values (these are all non-zero).
|
||||
#
|
||||
# It means that the ordinary least square model (with the optimal parameters) $\boldsymbol{\tilde{y}}$, corresponds to an orthogonal transformation of the output (or target) vector $\boldsymbol{y}$ by the vectors of the matrix $\boldsymbol{U}$.
|
||||
#
|
||||
#
|
||||
|
||||
# ## Further properties (important for our analyses later)
|
||||
#
|
||||
# Let us study again $\boldsymbol{X}^T\boldsymbol{X}$ in terms of our SVD,
|
||||
@@ -638,11 +598,9 @@ print(np.abs(C-B))
|
||||
# number of rows represents the number of data inputs. Note that in
|
||||
# other texts you may find the opposite notation. This has consequences
|
||||
# for the definition of for example the covariance matrix and its relation to the SVD.
|
||||
#
|
||||
#
|
||||
|
||||
# ## Meet the Covariance Matrix
|
||||
#
|
||||
#
|
||||
# Before we move on to a discussion of Ridge and Lasso regression, we want to show an important example of the above.
|
||||
#
|
||||
# We have already noted that the matrix $\boldsymbol{X}^T\boldsymbol{X}$ in ordinary
|
||||
@@ -666,8 +624,6 @@ print(np.abs(C-B))
|
||||
# the eigenvalues of the covariance matrix and the Hessian matrix in
|
||||
# terms of the singular values. Let us develop these arguments, as they will play an important role in our machine learning studies.
|
||||
#
|
||||
#
|
||||
#
|
||||
# Before we discuss the link between for example Ridge regression and the singular value decomposition, we need to remind ourselves about
|
||||
# the definition of the covariance and the correlation function. These are quantities that play a central role in machine learning methods.
|
||||
#
|
||||
@@ -709,7 +665,6 @@ print(np.abs(C-B))
|
||||
# **Scikit-Learn** or **nunmpy's** function calculate the covariance, this
|
||||
# quantity will be computed with a factor $1/(n-1)$.
|
||||
#
|
||||
#
|
||||
# The covariance takes values between zero and infinity and may thus
|
||||
# lead to problems with loss of numerical precision for particularly
|
||||
# large values. It is common to scale the covariance matrix by
|
||||
@@ -733,8 +688,6 @@ print(np.abs(C-B))
|
||||
|
||||
# In the above example this is the function we constructed using **pandas**.
|
||||
#
|
||||
#
|
||||
#
|
||||
# In our derivation of the various regression algorithms like **Ordinary Least Squares** or **Ridge regression**
|
||||
# we defined the design/feature matrix $\boldsymbol{X}$ as
|
||||
|
||||
@@ -865,8 +818,6 @@ print(C)
|
||||
#
|
||||
# The above procedure with **numpy** can be made more compact if we use **pandas**.
|
||||
#
|
||||
#
|
||||
#
|
||||
# We whow here how we can set up the correlation matrix using **pandas**, as done in this simple code
|
||||
|
||||
# In[7]:
|
||||
@@ -947,8 +898,6 @@ print(covariance_matrix)
|
||||
# drop these elements and construct a correlation
|
||||
# matrix without these elements.
|
||||
#
|
||||
#
|
||||
#
|
||||
# We can rewrite the covariance matrix in a more compact form in terms of the design/feature matrix $\boldsymbol{X}$ as
|
||||
|
||||
# $$
|
||||
@@ -983,12 +932,12 @@ print(covariance_matrix)
|
||||
# \end{bmatrix},
|
||||
# $$
|
||||
|
||||
# where we wrote $$\boldsymbol{C}[\boldsymbol{x}_0,\boldsymbol{x}_1] = \boldsymbol{C}[\boldsymbol{x}]$$ to indicate that this is the covariance of the vectors $\boldsymbol{x}$ of the design/feature matrix $\boldsymbol{X}$.
|
||||
# where we wrote $\boldsymbol{C}[\boldsymbol{x}_0,\boldsymbol{x}_1]=\boldsymbol{C}[\boldsymbol{x}]$ to indicate
|
||||
# that this is the covariance of the vectors $\boldsymbol{x}$ of the
|
||||
# design/feature matrix $\boldsymbol{X}$.
|
||||
#
|
||||
# It is easy to generalize this to a matrix $\boldsymbol{X}\in {\mathbb{R}}^{n\times p}$.
|
||||
#
|
||||
#
|
||||
#
|
||||
|
||||
# ## Linking with the SVD
|
||||
#
|
||||
# We saw earlier that
|
||||
@@ -1042,7 +991,6 @@ print(covariance_matrix)
|
||||
# $\boldsymbol{v}_i$ are hierarchically ordered by how much correlation they
|
||||
# encode from the columns of $\boldsymbol{X}$.
|
||||
#
|
||||
#
|
||||
# Note that these are also the eigenvectors and eigenvalues of the
|
||||
# Hessian matrix.
|
||||
#
|
||||
@@ -1060,7 +1008,6 @@ print(covariance_matrix)
|
||||
# self-adjoint, the singular values of $\boldsymbol{X}$ are equal to the
|
||||
# absolute value of the eigenvalues of $\boldsymbol{X}$.
|
||||
#
|
||||
#
|
||||
# For $\boldsymbol{X}\boldsymbol{X}^T$ we found
|
||||
|
||||
# $$
|
||||
@@ -1093,10 +1040,7 @@ print(covariance_matrix)
|
||||
# Since we will mainly be interested in the correlations among the features
|
||||
# of our data (the columns of $\boldsymbol{X}$, the quantity of interest for us are the non-zero singular
|
||||
# values and the column vectors of $\boldsymbol{V}$.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
|
||||
# ## Ridge and Lasso Regression
|
||||
#
|
||||
# Let us remind ourselves about the expression for the standard Mean Squared Error (MSE) which we used to define our cost function and the equations for the ordinary least squares (OLS) method, that is
|
||||
@@ -1184,15 +1128,12 @@ print(covariance_matrix)
|
||||
|
||||
# which can lead to singular matrices. However, with the SVD, we can always compute the inverse of the matrix $\boldsymbol{X}^T\boldsymbol{X}$.
|
||||
#
|
||||
#
|
||||
# We see that Ridge regression is nothing but the standard OLS with a
|
||||
# modified diagonal term added to $\boldsymbol{X}^T\boldsymbol{X}$. The consequences, in
|
||||
# particular for our discussion of the bias-variance tradeoff are rather
|
||||
# interesting. We will see that for specific values of $\lambda$, we may
|
||||
# even reduce the variance of the optimal parameters $\boldsymbol{\beta}$. These topics and other related ones, will be discussed after the more linear algebra oriented analysis here.
|
||||
#
|
||||
#
|
||||
#
|
||||
# Using our insights about the SVD of the design matrix $\boldsymbol{X}$
|
||||
# We have already analyzed the OLS solutions in terms of the eigenvectors (the columns) of the right singular value matrix $\boldsymbol{U}$ as
|
||||
|
||||
@@ -1208,8 +1149,6 @@ print(covariance_matrix)
|
||||
|
||||
# with the vectors $\boldsymbol{u}_j$ being the columns of $\boldsymbol{U}$ from the SVD of the matrix $\boldsymbol{X}$.
|
||||
#
|
||||
#
|
||||
#
|
||||
# Since $\lambda \geq 0$, it means that compared to OLS, we have
|
||||
|
||||
# $$
|
||||
@@ -1224,8 +1163,6 @@ print(covariance_matrix)
|
||||
#
|
||||
# For small eigenvalues $\sigma_i$ it means that their contributions become less important, a fact which can be used to reduce the number of degrees of freedom. More about this when we have covered the material on a statistical interpretation of various linear regression methods.
|
||||
#
|
||||
#
|
||||
#
|
||||
# For the sake of simplicity, let us assume that the design matrix is orthonormal, that is
|
||||
|
||||
# $$
|
||||
@@ -1250,8 +1187,6 @@ print(covariance_matrix)
|
||||
#
|
||||
# We will come back to more interpreations after we have gone through some of the statistical analysis part.
|
||||
#
|
||||
#
|
||||
#
|
||||
# Using the matrix-vector expression for Lasso regression and dropping the parameter $1/n$ in front of the standard mean squared error equation, we have the following **cost** function
|
||||
|
||||
# $$
|
||||
@@ -1278,10 +1213,6 @@ print(covariance_matrix)
|
||||
|
||||
# This equation does not lead to a nice analytical equation as in Ridge regression or ordinary least squares. This equation can however be solved by using standard convex optimization algorithms using for example the Python package [CVXOPT](https://cvxopt.org/). We will discuss this later.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# Let us assume that our design matrix is given by unit (identity) matrix, that is a square diagonal matrix with ones only along the
|
||||
# diagonal. In this case we have an equal number of rows and columns $n=p$.
|
||||
#
|
||||
@@ -1331,7 +1262,6 @@ print(covariance_matrix)
|
||||
|
||||
# Plotting these results ([figure in handwritten notes for week 36](https://github.com/CompPhysics/MachineLearning/blob/master/doc/HandWrittenNotes/2021/NotesSeptember9.pdf)) shows clearly that Lasso regression suppresses (sets to zero) values of $\beta_i$ for specific values of $\lambda$. Ridge regression reduces on the other hand the values of $\beta_i$ as function of $\lambda$.
|
||||
#
|
||||
#
|
||||
# As another examples,
|
||||
# let us assume we have a data set with outputs/targets given by the vector
|
||||
|
||||
@@ -1347,7 +1277,6 @@ print(covariance_matrix)
|
||||
|
||||
# meaning that we have two features and two unknown parameters $\beta_0$ and $\beta_1$ to be determined either by ordinary least squares, Ridge or Lasso regression.
|
||||
#
|
||||
#
|
||||
# For ordinary least squares (OLS) we know that the optimal solution is
|
||||
|
||||
# $$
|
||||
@@ -1362,7 +1291,6 @@ print(covariance_matrix)
|
||||
|
||||
# The code which implements this simpler case is presented after the discussion of Ridge and Lasso.
|
||||
#
|
||||
#
|
||||
# For Ridge regression we have
|
||||
|
||||
# $$
|
||||
@@ -1380,8 +1308,6 @@ print(covariance_matrix)
|
||||
#
|
||||
# To see this, let us write the cost function for Ridge regression.
|
||||
#
|
||||
#
|
||||
#
|
||||
# We define the MSE without the $1/n$ factor and have then, using that
|
||||
|
||||
# $$
|
||||
@@ -1412,7 +1338,6 @@ print(covariance_matrix)
|
||||
|
||||
# which gives $\lambda=4.571$ and $\beta_0=0.933$ and $\beta_1=0.359$.
|
||||
#
|
||||
#
|
||||
# For Lasso we need now, keeping a constraint on $\vert\beta_0\vert+\vert\beta_1\vert=1$, to take the derivative of the absolute values of $\beta_0$
|
||||
# and $\beta_1$. This gives us the following derivatives of the cost function
|
||||
|
||||
@@ -1465,7 +1390,6 @@ print(covariance_matrix)
|
||||
|
||||
# Using the constraint on $\beta_0$ and $\beta_1$ we can then find the optimal value of $\lambda$ for the different cases. We leave this as an exercise to you.
|
||||
#
|
||||
#
|
||||
# Here we set up the OLS, Ridge and Lasso functionality in order to study the above example. Note that here we have opted for a set of values of $\lambda$, meaning that we need to perform a search in order to find the optimal values.
|
||||
#
|
||||
# First we study and compare the OLS and Ridge results. The next code compares all three methods.
|
||||
@@ -1694,285 +1618,8 @@ plt.show()
|
||||
# for the optimal values of $\lambda$. We will see this throughout these
|
||||
# series of lectures.
|
||||
#
|
||||
#
|
||||
# As a small addendum, we note that you can also solve this problem using the convex optimization package [CVXOPT](https://cvxopt.org/examples/mlbook/l1regls.html). This requires, in addition to having installed **CVXOPT**, you need to download the file *l1regl.py*.
|
||||
# The following code example solves the simpler problem we discussed above, where we have added the latter python file.
|
||||
|
||||
# In[12]:
|
||||
|
||||
|
||||
from cvxopt import matrix, spdiag, mul, div, sqrt, normal, setseed
|
||||
from cvxopt import blas, lapack, solvers, sparse, spmatrix
|
||||
import math
|
||||
|
||||
try:
|
||||
import mosek
|
||||
import sys
|
||||
__MOSEK = True
|
||||
except: __MOSEK = False
|
||||
|
||||
if __MOSEK:
|
||||
|
||||
def l1regls_mosek(A, b):
|
||||
"""
|
||||
|
||||
Returns the solution of l1-norm regularized least-squares problem
|
||||
|
||||
minimize || A*x - b ||_2^2 + e'*u
|
||||
|
||||
subject to -u <= x <= u
|
||||
|
||||
"""
|
||||
|
||||
m, n = A.size
|
||||
|
||||
env = mosek.Env()
|
||||
task = env.Task(0,0)
|
||||
task.set_Stream(mosek.streamtype.log, lambda x: sys.stdout.write(x))
|
||||
|
||||
task.appendvars( 2*n) # number of variables
|
||||
task.appendcons( 2*n) # number of constraints
|
||||
|
||||
# input quadratic objective
|
||||
Q = matrix(0.0, (n,n))
|
||||
blas.syrk(A, Q, alpha = 2.0, trans='T')
|
||||
|
||||
I = []
|
||||
for i in range(n):
|
||||
I.extend(range(i,n))
|
||||
|
||||
J = []
|
||||
for i in range(n):
|
||||
J.extend((n-i)*[i])
|
||||
|
||||
task.putqobj(I, J, list(Q[matrix(I) + matrix(J)*n]))
|
||||
task.putclist(range(2*n), list(-2*A.T*b) + n*[1.0]) # setup linear objective
|
||||
|
||||
# input constraint matrix row by row
|
||||
for i in range(n):
|
||||
task.putarow( i, [i, n+i], [1.0, -1.0])
|
||||
task.putarow( n+i, [i, n+i], [1.0, 1.0])
|
||||
|
||||
# setup bounds on constraints
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
0, n, n*[mosek.boundkey.up], n*[0.0], n*[0.0])
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
n, 2*n, n*[mosek.boundkey.lo], n*[0.0], n*[0.0])
|
||||
|
||||
# setup variable bounds
|
||||
task.putboundslice(mosek.accmode.var,
|
||||
0, 2*n, 2*n*[mosek.boundkey.fr], 2*n*[0.0], 2*n*[0.0])
|
||||
|
||||
# optimize the task
|
||||
task.putobjsense(mosek.objsense.minimize)
|
||||
task.optimize()
|
||||
task.solutionsummary(mosek.streamtype.log)
|
||||
x = n*[0.0]
|
||||
task.getsolutionslice(mosek.soltype.itr, mosek.solitem.xx, 0, n, x)
|
||||
|
||||
return matrix(x)
|
||||
|
||||
def l1regls_mosek2(A, b):
|
||||
"""
|
||||
|
||||
Returns the solution of l1-norm regularized least-squares problem
|
||||
|
||||
minimize w'*w + e'*u
|
||||
|
||||
subject to -u <= x <= u
|
||||
|
||||
A*x - w = b
|
||||
|
||||
"""
|
||||
|
||||
m, n = A.size
|
||||
|
||||
env = mosek.Env()
|
||||
task = env.Task(0,0)
|
||||
task.set_Stream(mosek.streamtype.log, lambda x: sys.stdout.write(x))
|
||||
|
||||
task.appendvars(2*n + m) # number of variables
|
||||
task.appendcons(2*n + m) # number of constraints
|
||||
|
||||
# input quadratic objective
|
||||
task.putqobj(range(2*n,2*n+m), range(2*n,2*n+m), m*[2.0])
|
||||
|
||||
task.putclist(range(2*n+m), n*[0.0] + n*[1.0] + m*[0.0]) # setup linear objective
|
||||
|
||||
# input constraint matrix row by row
|
||||
for i in range(n):
|
||||
task.putarow( i, [i, n+i], [1.0, -1.0])
|
||||
task.putarow( n+i, [i, n+i], [1.0, 1.0])
|
||||
|
||||
for i in range(m):
|
||||
task.putarow( 2*n+i, range(n) + [2*n+i], list(A[i,:]) + [-1.0])
|
||||
|
||||
# setup bounds on constraints
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
0, n, n*[mosek.boundkey.up], n*[0.0], n*[0.0])
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
n, 2*n, n*[mosek.boundkey.lo], n*[0.0], n*[0.0])
|
||||
task.putboundslice(mosek.accmode.con,
|
||||
2*n, 2*n+m, m*[mosek.boundkey.fx], list(b), list(b))
|
||||
|
||||
# setup variable bounds
|
||||
task.putboundslice(mosek.accmode.var, 0, 2*n+m, (2*n+m)*[mosek.boundkey.fr],
|
||||
(2*n+m)*[0.0], (2*n+m)*[0.0])
|
||||
|
||||
# optimize the task
|
||||
task.putobjsense(mosek.objsense.minimize)
|
||||
task.optimize()
|
||||
task.solutionsummary(mosek.streamtype.log)
|
||||
x = n*[0.0]
|
||||
task.getsolutionslice(mosek.soltype.itr, mosek.solitem.xx, 0, n, x)
|
||||
|
||||
return matrix(x)
|
||||
|
||||
def l1regls(A, b):
|
||||
"""
|
||||
|
||||
Returns the solution of l1-norm regularized least-squares problem
|
||||
|
||||
minimize || A*x - b ||_2^2 + || x ||_1.
|
||||
|
||||
"""
|
||||
|
||||
m, n = A.size
|
||||
q = matrix(1.0, (2*n,1))
|
||||
q[:n] = -2.0 * A.T * b
|
||||
|
||||
def P(u, v, alpha = 1.0, beta = 0.0 ):
|
||||
"""
|
||||
v := alpha * 2.0 * [ A'*A, 0; 0, 0 ] * u + beta * v
|
||||
"""
|
||||
v *= beta
|
||||
v[:n] += alpha * 2.0 * A.T * (A * u[:n])
|
||||
|
||||
|
||||
def G(u, v, alpha=1.0, beta=0.0, trans='N'):
|
||||
"""
|
||||
v := alpha*[I, -I; -I, -I] * u + beta * v (trans = 'N' or 'T')
|
||||
"""
|
||||
|
||||
v *= beta
|
||||
v[:n] += alpha*(u[:n] - u[n:])
|
||||
v[n:] += alpha*(-u[:n] - u[n:])
|
||||
|
||||
h = matrix(0.0, (2*n,1))
|
||||
|
||||
|
||||
# Customized solver for the KKT system
|
||||
#
|
||||
# [ 2.0*A'*A 0 I -I ] [x[:n] ] [bx[:n] ]
|
||||
# [ 0 0 -I -I ] [x[n:] ] = [bx[n:] ].
|
||||
# [ I -I -D1^-1 0 ] [zl[:n]] [bzl[:n]]
|
||||
# [ -I -I 0 -D2^-1 ] [zl[n:]] [bzl[n:]]
|
||||
#
|
||||
# where D1 = W['di'][:n]**2, D2 = W['di'][:n]**2.
|
||||
#
|
||||
# We first eliminate zl and x[n:]:
|
||||
#
|
||||
# ( 2*A'*A + 4*D1*D2*(D1+D2)^-1 ) * x[:n] =
|
||||
# bx[:n] - (D2-D1)*(D1+D2)^-1 * bx[n:] +
|
||||
# D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] -
|
||||
# D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:]
|
||||
#
|
||||
# x[n:] = (D1+D2)^-1 * ( bx[n:] - D1*bzl[:n] - D2*bzl[n:] )
|
||||
# - (D2-D1)*(D1+D2)^-1 * x[:n]
|
||||
#
|
||||
# zl[:n] = D1 * ( x[:n] - x[n:] - bzl[:n] )
|
||||
# zl[n:] = D2 * (-x[:n] - x[n:] - bzl[n:] ).
|
||||
#
|
||||
# The first equation has the form
|
||||
#
|
||||
# (A'*A + D)*x[:n] = rhs
|
||||
#
|
||||
# and is equivalent to
|
||||
#
|
||||
# [ D A' ] [ x:n] ] = [ rhs ]
|
||||
# [ A -I ] [ v ] [ 0 ].
|
||||
#
|
||||
# It can be solved as
|
||||
#
|
||||
# ( A*D^-1*A' + I ) * v = A * D^-1 * rhs
|
||||
# x[:n] = D^-1 * ( rhs - A'*v ).
|
||||
|
||||
S = matrix(0.0, (m,m))
|
||||
Asc = matrix(0.0, (m,n))
|
||||
v = matrix(0.0, (m,1))
|
||||
|
||||
def Fkkt(W):
|
||||
|
||||
# Factor
|
||||
#
|
||||
# S = A*D^-1*A' + I
|
||||
#
|
||||
# where D = 2*D1*D2*(D1+D2)^-1, D1 = d[:n]**-2, D2 = d[n:]**-2.
|
||||
|
||||
d1, d2 = W['di'][:n]**2, W['di'][n:]**2
|
||||
|
||||
# ds is square root of diagonal of D
|
||||
ds = math.sqrt(2.0) * div( mul( W['di'][:n], W['di'][n:]),
|
||||
sqrt(d1+d2) )
|
||||
d3 = div(d2 - d1, d1 + d2)
|
||||
|
||||
# Asc = A*diag(d)^-1/2
|
||||
Asc = A * spdiag(ds**-1)
|
||||
|
||||
# S = I + A * D^-1 * A'
|
||||
blas.syrk(Asc, S)
|
||||
S[::m+1] += 1.0
|
||||
lapack.potrf(S)
|
||||
|
||||
def g(x, y, z):
|
||||
|
||||
x[:n] = 0.5 * ( x[:n] - mul(d3, x[n:]) +
|
||||
mul(d1, z[:n] + mul(d3, z[:n])) - mul(d2, z[n:] -
|
||||
mul(d3, z[n:])) )
|
||||
x[:n] = div( x[:n], ds)
|
||||
|
||||
# Solve
|
||||
#
|
||||
# S * v = 0.5 * A * D^-1 * ( bx[:n] -
|
||||
# (D2-D1)*(D1+D2)^-1 * bx[n:] +
|
||||
# D1 * ( I + (D2-D1)*(D1+D2)^-1 ) * bzl[:n] -
|
||||
# D2 * ( I - (D2-D1)*(D1+D2)^-1 ) * bzl[n:] )
|
||||
|
||||
blas.gemv(Asc, x, v)
|
||||
lapack.potrs(S, v)
|
||||
|
||||
# x[:n] = D^-1 * ( rhs - A'*v ).
|
||||
blas.gemv(Asc, v, x, alpha=-1.0, beta=1.0, trans='T')
|
||||
x[:n] = div(x[:n], ds)
|
||||
|
||||
# x[n:] = (D1+D2)^-1 * ( bx[n:] - D1*bzl[:n] - D2*bzl[n:] )
|
||||
# - (D2-D1)*(D1+D2)^-1 * x[:n]
|
||||
x[n:] = div( x[n:] - mul(d1, z[:n]) - mul(d2, z[n:]), d1+d2 ) - mul( d3, x[:n] )
|
||||
|
||||
# zl[:n] = D1^1/2 * ( x[:n] - x[n:] - bzl[:n] )
|
||||
# zl[n:] = D2^1/2 * ( -x[:n] - x[n:] - bzl[n:] ).
|
||||
z[:n] = mul( W['di'][:n], x[:n] - x[n:] - z[:n] )
|
||||
z[n:] = mul( W['di'][n:], -x[:n] - x[n:] - z[n:] )
|
||||
|
||||
return g
|
||||
|
||||
return solvers.coneqp(P, q, G, h, kktsolver = Fkkt)['x'][:n]
|
||||
|
||||
|
||||
# Then we call the above functions and solve the problem, as done here
|
||||
|
||||
# In[ ]:
|
||||
|
||||
|
||||
from cvxopt import matrix, normal
|
||||
|
||||
X = matrix( [ [ 2, 0, 1], [0, 1, 3]])
|
||||
y = matrix( [4, 2, 3])
|
||||
x = l1regls(X,y)
|
||||
|
||||
|
||||
# **More text will be added to this example.**
|
||||
#
|
||||
# ## Linking the regression analysis with a statistical interpretation
|
||||
#
|
||||
# We will now couple the discussions of ordinary least squares, Ridge
|
||||
@@ -1983,7 +1630,6 @@ x = l1regls(X,y)
|
||||
# parameter can reduce considerably the variance of the parameters
|
||||
# $\beta$.
|
||||
#
|
||||
#
|
||||
# The
|
||||
# advantage of doing linear regression is that we actually end up with
|
||||
# analytical expressions for several statistical quantities.
|
||||
@@ -1991,7 +1637,6 @@ x = l1regls(X,y)
|
||||
# derive quantities like the variance and other expectation values in a
|
||||
# rather straightforward way.
|
||||
#
|
||||
#
|
||||
# It is assumed that $\varepsilon_i
|
||||
# \sim \mathcal{N}(0, \sigma^2)$ and the $\varepsilon_{i}$ are
|
||||
# independent, i.e.:
|
||||
@@ -2015,8 +1660,6 @@ x = l1regls(X,y)
|
||||
# notation above $\mathbf{X}_{i,\ast}$ means that we are looking at the
|
||||
# row number $i$ and perform a sum over all values $p$.
|
||||
#
|
||||
#
|
||||
#
|
||||
# The assumption we have made here can be summarized as (and this is going to be useful when we discuss the bias-variance trade off)
|
||||
# that there exists a function $f(\boldsymbol{x})$ and a normal distributed error $\boldsymbol{\varepsilon}\sim \mathcal{N}(0, \sigma^2)$
|
||||
# which describe our data
|
||||
@@ -2063,7 +1706,6 @@ x = l1regls(X,y)
|
||||
# Hence, $y_i \sim \mathcal{N}( \mathbf{X}_{i, \ast} \, \boldsymbol{\beta}, \sigma^2)$, that is $\boldsymbol{y}$ follows a normal distribution with
|
||||
# mean value $\boldsymbol{X}\boldsymbol{\beta}$ and variance $\sigma^2$ (not be confused with the singular values of the SVD).
|
||||
#
|
||||
#
|
||||
# With the OLS expressions for the parameters $\boldsymbol{\beta}$ we can evaluate the expectation value
|
||||
|
||||
# $$
|
||||
@@ -2107,7 +1749,6 @@ x = l1regls(X,y)
|
||||
# $\boldsymbol{\sigma}^2 (\boldsymbol{\beta}_j ) = \boldsymbol{\sigma}^2 [(\mathbf{X}^{T} \mathbf{X})^{-1}]_{jj} $. This may be used to
|
||||
# construct a confidence interval for the estimates.
|
||||
#
|
||||
#
|
||||
# In a similar way, we can obtain analytical expressions for say the
|
||||
# expectation values of the parameters $\boldsymbol{\beta}$ and their variance
|
||||
# when we employ Ridge regression, allowing us again to define a confidence interval.
|
||||
@@ -2137,10 +1778,8 @@ x = l1regls(X,y)
|
||||
|
||||
# The difference is non-negative definite since each component of the
|
||||
# matrix product is non-negative definite.
|
||||
# This means the variance we obtain with the standard OLS will always for $\lambda > 0$ be larger than the variance of $\boldsymbol{\beta}$ obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below.
|
||||
#
|
||||
#
|
||||
#
|
||||
# This means the variance we obtain with the standard OLS will always for $\lambda > 0$ be larger than the variance of $\boldsymbol{\beta}$ obtained with the Ridge estimator. This has interesting consequences when we discuss the so-called bias-variance trade-off below.
|
||||
|
||||
# ## Deriving OLS from a probability distribution
|
||||
#
|
||||
# Our basic assumption when we derived the OLS equations was to assume
|
||||
@@ -2193,21 +1832,18 @@ x = l1regls(X,y)
|
||||
# likelihood of a domain of events $\boldsymbol{D}$ given a set of parameters
|
||||
# $\boldsymbol{\beta}$.
|
||||
#
|
||||
#
|
||||
# In statistics, maximum likelihood estimation (MLE) is a method of
|
||||
# estimating the parameters of an assumed probability distribution,
|
||||
# given some observed data. This is achieved by maximizing a likelihood
|
||||
# function so that, under the assumed statistical model, the observed
|
||||
# data is the most probable.
|
||||
#
|
||||
#
|
||||
# We will assume here that our events are given by the above Gaussian
|
||||
# distribution and we will determine the optimal parameters $\beta$ by
|
||||
# maximizing the above PDF. However, computing the derivatives of a
|
||||
# product function is cumbersome and can easily lead to overflow and/or
|
||||
# underflowproblems, with potentials for loss of numerical precision.
|
||||
#
|
||||
#
|
||||
# In practice, it is more convenient to maximize the logarithm of the
|
||||
# PDF because it is a monotonically increasing function of the argument.
|
||||
# Alternatively, and this will be our option, we will minimize the
|
||||
@@ -2217,9 +1853,6 @@ x = l1regls(X,y)
|
||||
# Note also that maximization/minimization of the logarithm of the PDF
|
||||
# is equivalent to the maximization/minimization of the function itself.
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# We could now define a new cost function to minimize, namely the negative logarithm of the above PDF
|
||||
|
||||
# $$
|
||||
@@ -2246,7 +1879,6 @@ x = l1regls(X,y)
|
||||
|
||||
# Before we make a similar analysis for Ridge and Lasso regression, we need a short reminder on statistics.
|
||||
#
|
||||
#
|
||||
# A central theorem in statistics is Bayes' theorem. This theorem plays a similar role as the good old Pythagoras' theorem in geometry.
|
||||
# Bayes' theorem is extremely simple to derive. But to do so we need some basic axioms from statistics.
|
||||
#
|
||||
@@ -2271,8 +1903,6 @@ x = l1regls(X,y)
|
||||
#
|
||||
# If we have independent events then $p(X,Y)=p(X)p(Y)$.
|
||||
#
|
||||
#
|
||||
#
|
||||
# The marginal probability is defined in terms of only one of the set of variables $X,Y$. For a discrete probability we have
|
||||
|
||||
# $$
|
||||
@@ -2299,7 +1929,6 @@ x = l1regls(X,y)
|
||||
|
||||
# which is Bayes' theorem. It allows us to evaluate the uncertainty in in $X$ after we have observed $Y$. We can easily interchange $X$ with $Y$.
|
||||
#
|
||||
#
|
||||
# The quantity $p(Y\vert X)$ on the right-hand side of the theorem is
|
||||
# evaluated for the observed data $Y$ and can be viewed as a function of
|
||||
# the parameter space represented by $X$. This function is not
|
||||
@@ -2312,7 +1941,6 @@ x = l1regls(X,y)
|
||||
#
|
||||
# Let us try to illustrate Bayes' theorem through an example.
|
||||
#
|
||||
#
|
||||
# Let us suppose that you are undergoing a series of mammography scans
|
||||
# in order to rule out possible breast cancer cases. We define the
|
||||
# sensitivity for a positive event by the variable $X$. It takes binary
|
||||
@@ -2342,8 +1970,6 @@ x = l1regls(X,y)
|
||||
|
||||
# instead of $p(X=1\vert Y=1)$.
|
||||
#
|
||||
#
|
||||
#
|
||||
# If we look at various national surveys on breast cancer, the general
|
||||
# likelihood of developing breast cancer is a very small number. Let us
|
||||
# assume that the prior probability in the population as a whole is
|
||||
@@ -2384,9 +2010,7 @@ x = l1regls(X,y)
|
||||
# $$
|
||||
|
||||
# That is, in case of a positive test, there is only a $3\%$ chance of having breast cancer!
|
||||
#
|
||||
#
|
||||
#
|
||||
|
||||
# ## Bayes' Theorem and Ridge and Lasso Regression
|
||||
#
|
||||
# Hitherto we have discussed Ridge and Lasso regression in terms of a
|
||||
@@ -2396,7 +2020,6 @@ x = l1regls(X,y)
|
||||
#
|
||||
# Before we proceed let us perform a Ridge, Lasso and OLS analysis of a polynomial fit.
|
||||
#
|
||||
#
|
||||
# We will play around with a study of the values for the optimal
|
||||
# parameters $\boldsymbol{\beta}$ using OLS, Ridge and Lasso regression. For
|
||||
# OLS, you will notice as function of the noise and polynomial degree,
|
||||
@@ -2408,7 +2031,7 @@ x = l1regls(X,y)
|
||||
# typically be reduced, providing thereby less fluctuations from one
|
||||
# order to another one.
|
||||
|
||||
# In[ ]:
|
||||
# In[12]:
|
||||
|
||||
|
||||
import numpy as np
|
||||
@@ -2497,7 +2120,7 @@ plt.show()
|
||||
# quench the fluctuations in the parameters of $\beta_i$ which have a
|
||||
# large variance (normally for higher orders in the polynomial).
|
||||
|
||||
# In[ ]:
|
||||
# In[13]:
|
||||
|
||||
|
||||
import numpy as np
|
||||
@@ -2546,11 +2169,8 @@ for i in range(nlambdas):
|
||||
# noise. Here we recommend to use $\sigma^2=1$ as variance for the
|
||||
# added noise (which follows a normal distribution with mean value zero).
|
||||
# Comment your results. If you have a large noise term, do the parameters $\beta_j$ vary more as function
|
||||
# model complexity? And what about their variance?
|
||||
#
|
||||
#
|
||||
#
|
||||
#
|
||||
# model complexity? And what about their variance?
|
||||
|
||||
# ## Linking Bayes' Theorem with Ridge and Lasso Regression
|
||||
#
|
||||
# We have seen that Ridge regression suppresses those features which
|
||||
@@ -2585,8 +2205,6 @@ for i in range(nlambdas):
|
||||
|
||||
# We have a model for $p(\boldsymbol{D}\vert\boldsymbol{\beta})$ but need one for the **prior** $p(\boldsymbol{\beta}$!
|
||||
#
|
||||
#
|
||||
#
|
||||
# With the posterior probability defined by a likelihood which we have
|
||||
# already modeled and an unknown prior, we are now ready to make
|
||||
# additional models for the prior.
|
||||
@@ -2594,7 +2212,6 @@ for i in range(nlambdas):
|
||||
# We can, based on our discussions of the variance of $\boldsymbol{\beta}$ and
|
||||
# the mean value, assume that the prior for the values $\boldsymbol{\beta}$ is
|
||||
# given by a Gaussian with mean value zero and variance $\tau^2$, that
|
||||
# is
|
||||
|
||||
# $$
|
||||
# p(\boldsymbol{\beta})=\prod_{j=0}^{p-1}\exp{\left(-\frac{\beta_j^2}{2\tau^2}\right)}.
|
||||
@@ -2623,7 +2240,6 @@ for i in range(nlambdas):
|
||||
|
||||
# which is our Ridge cost function! Nice, isn't it?
|
||||
#
|
||||
#
|
||||
# To derive the Lasso cost function, we simply replace the Gaussian prior with an exponential distribution ([Laplace in this case](https://en.wikipedia.org/wiki/Laplace_distribution)) with zero mean value, that is
|
||||
|
||||
# $$
|
||||
@@ -2652,7 +2268,6 @@ for i in range(nlambdas):
|
||||
|
||||
# which is our Lasso cost function!
|
||||
#
|
||||
#
|
||||
# Plotting these prior functions shows us that we can use the parameter
|
||||
# $\lambda$ to shrink or increase the role of a given parameter
|
||||
# $\beta_j$. The variance for the Laplace distribution is
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 9.1 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 11 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 14 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 10 KiB |
+1477
-787
File diff suppressed because it is too large
Load Diff
Reference in New Issue
Block a user