In [3]:
using LinearAlgebra, Plots, Random, SparseArrays
In [ ]:
In [2]:
"""
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 = qr(B')
    R = qrf.R
    @assert all(diag(R) != 0) "B is not full-rank"
    
    Q = qrf.Q*Matrix(I,n,n) # this forms the full matrix Q. 
    Q1 = view(Q,:,1:p)     # view extracts a submatrix
    Q2 = view(Q,:,(p+1):n) 
    
    # form a solution
    xc = Q1*(R'\d)
    # solve the new LS problem
    y = (A*Q2)\(b-A*(xc))
    # recombine the solution
    x =xc + Q2*y
    return x
end
Out[2]:
conls
In [ ]:
In [10]:
teams = ["duke","miami","unc","uva","vt"]
data = [ # team 1 team 2, team 1 pts, team 2 pts
    1 2  7 52 # duke played Miami and lost 7 to 52 
    1 3 21 24 # duke played unc and lost 21 to 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]
    Pi = data[g,3]
    Pj = data[g,4]
  
    G[g,i] = 1
    G[g,j] = -1
    p[g] = Pi - Pj
end
s1 = conls(G,p,ones(1,nteams),[200.0])
Out[10]:
5×1 Matrix{Float64}:
 15.199999999999992
 58.2
 31.999999999999993
 36.599999999999994
 58.0
In [5]:
sum(s1)
Out[5]:
99.99999999999999
In [6]:
G
Out[6]:
10×5 Matrix{Float64}:
 1.0  -1.0   0.0   0.0   0.0
 1.0   0.0  -1.0   0.0   0.0
 1.0   0.0   0.0  -1.0   0.0
 1.0   0.0   0.0   0.0  -1.0
 0.0   1.0  -1.0   0.0   0.0
 0.0   1.0   0.0  -1.0   0.0
 0.0   1.0   0.0   0.0  -1.0
 0.0   0.0   1.0  -1.0   0.0
 0.0   0.0   1.0   0.0  -1.0
 0.0   0.0   0.0   1.0  -1.0
In [7]:
for i=1:nteams
    println(teams[i], " ", s1[i])
end

duke -4.800000000000004
miami 38.2
unc 11.999999999999996
uva 16.599999999999998
vt 38.0
In [18]:
t1 = zeros(1,nteams)
t1[1] = 1
s3 = conls(G,p,t1,[Float64(pi)])
Out[18]:
5×1 Matrix{Float64}:
  3.141592653589793
 46.1415926535898
 19.9415926535898
 24.541592653589802
 45.9415926535898

Look at the residual

In [ ]:
G*s1 - p
In [ ]:
In [ ]:
sum(s1)
In [11]:
"""
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[11]:
conls_lagrange
In [14]:
s2,lam = conls_lagrange(G,p,ones(1,nteams),[200.0])
Out[14]:
([15.2, 58.2, 31.999999999999996, 36.6, 58.0], [0.0])
In [15]:
s1 - s2
Out[15]:
5×1 Matrix{Float64}:
 -7.105427357601002e-15
  0.0
 -3.552713678800501e-15
 -7.105427357601002e-15
  0.0
In [ ]: