In [1]:
"""
conls
=====
Solve a constrained least squares problem
conls(A,b,B,d) -> x
solves
minimize ||b - Ax||_2
subject to Bx = d
when B is p-by-n, p < n, and rank(B) = p, by using the QR factorization
of B' to find a basis for the nullspace of B.
"""
function conls(A,b,B,d)
m,n = size(A)
p = size(B,1)
@assert size(B,2) == n "incompatible dimensions for A and B"
# Step 1, QR of B' to get the nullspace
qrf = qrfact(B')
R = qrf[:R]
@assert all(diag(R) != 0) "B is not full-rank"
Q = full(qrf[:Q], thin=false)
Q1 = view(Q,:,1:p) # view extracts a submatrix
Q2 = view(Q,:,(p+1):n)
# form a solution
y = R'\d
# solve the new LS problem
z = (A*Q2)\(b-A*(Q1*y))
# recombine the solution
x = Q1*y + Q2*z
return x
end
Out[1]:
In [2]:
teams = ["duke","miami","unc","uva","vt"]
data = [ # team 1 team 2, team 1 pts, team 2 pts
1 2 7 52
1 3 21 24
1 4 7 38
1 5 0 45
2 3 34 16
2 4 25 17
2 5 27 7
3 4 7 5
3 5 3 30
4 5 14 52
]
ngames = size(data,1)
nteams = length(teams)
G = zeros(ngames, nteams)
p = zeros(ngames, 1)
for g=1:ngames
i = data[g,1]
j = data[g,2]
Si = data[g,3]
Sj = data[g,4]
G[g,i] = 1
G[g,j] = -1
p[g] = Si - Sj
end
s1 = conls(G,p,ones(1,nteams),[1.0])
Out[2]:
In [3]:
sum(s1)
Out[3]:
In [4]:
"""
conls_lagrange
=====
Solve a constrained least squares problem
conls(A,b,B,d) -> x, lam
solves
minimize ||b - Ax||_2
subject to Bx = d
when B is p-by-n, p < n, and rank(B) = p, by using the augmented
system of normal equations
[ A'*A B' ] [ x ] = [ A'*b]
[ B 0 ] [ lam ] = [ d ]
and returns both x and the lagrange multipliers
"""
function conls_lagrange(A,b,B,d)
m,n = size(A)
p = size(B,1)
@assert size(B,2) == n "incompatible dimensions for A and B"
# Step 1, form the block system
M = [A'*A B'; B zeros(p,p) ]
# Step 2, solve
z = M\[A'*b; d]
# Step 3, extract
x = z[1:n]
lam = z[n+1:end]
return x, lam
end
Out[4]:
In [5]:
s2,lam = conls_lagrange(G,p,ones(1,nteams),[1.0])
Out[5]: