Large Scale Optimization for Data Science, Jacek Gondzio 1 Matlab Assignment 1 (MA1)
This Assignment will be assessed and its result will count into the overall mark in LSO for DS. The following marking scheme will be used:
Problem1: 4+5+5+6=20marks; Problem2: 5+5+5=15marks;
Problem 3: 20 marks; Problem 4: 30 + 15 = 45 marks.
Maximum achievable: 100 marks.
MA1 is due before 2pm on Thursday 11th February.
Please download all the folders and the MATLAB routines contained in MA1, complete the missing parts in the routines and submit your solutions via Learn. (Please gather all your folders in one single folder.)
Please read the comments written in the routines for more details on the formats and types of INPUT and OUTPUT.
For each problem use the file Prob-xx-.m (and instructions therein) to check the MATLAB routines you produced.
Problem 1
Consider the following quadratic optimization problem
min x1 +2×2 + 1×21 + 3×2
s.t.
It has an optimal solution xˆ1 = 1 and xˆ2 = 0.
(i) Implement the routine
22 x1 + x2 = 1
x1 ≥0,×2 ≥0.
function [A,b,c,Q] = getInput()
to define the matrices A, Q and the vectors b, c when the optimization problem is written in the form of the standard quadratic program
min cTx+1xTQx 2
s.t. Ax = b x ≥ 0.
(1)
(ii) Write down the first order optimality conditions (FOC) for the associated barrier problem and implement the routine which computes the violation (the right-hand-side) of FOC:
function [xi p, xi d, xi mu] = getRhs(A, b, c, Q, mu, x, y, s)
(iii) Implement the routine
function [ dx, dy, ds ]= newtonDirections( A, Q, xi p, xi d, xi mu, x, s)
which determines the Newton directions.
(iv) Now, you should be able to find primal-dual strictly feasible solutions for the first order optimality conditions for any value of the barrier parameter μ. To do this, use the routines above and implement the Newton method in the routine
function x = feasibleSolution(mu,x0,y0,s0) which computes (x, y, s) such that
∥getRhs(A,b,c,Q,mu,x,y,s)∥2 <10−3
for a fixed INPUT μ and starting point (x0,y0,s0).
Observation: Do the solutions x = (x1(μ),x2(μ)) seem to converge to xˆ as μ is reduced?
Large Scale Optimization for Data Science, Jacek Gondzio 2
Problem 2
Recall the Sherman-Morrison-Woodbury formula. Let U, V ∈ Rn×k and define a rank-k update of matrix A ∈ Rn×n
A ̄ = A + U V T .
A ̄−1 =(A+UVT)−1 =A−1 −A−1U(Ik +VTA−1U)−1VTA−1.
(i) Implement the technique in the routine
function x = SMW(A, U, V, b)
which takes as an INPUT:
A∈Rn×n,U∈Rn×k,V ∈Rn×k andb∈Rn×r (r=1foracolumnvector) and produces an OUTPUT:
x, a solution of equation (A + U V T ) x = b.
(HINT: you do not need to compute the inverse of A but just to solve linear systems involving A.)
(ii) Modify the MATLAB functionSMW2.mwhich usesSMWto compute the inverse of A+UVT (HINT:
compute an inverse solving linear systems), where
251
161
152
A ̄= 2 1 1=I6+UVT,
1 2 1 112
UT=1 1 1 0 0 0andVT=1 5 1 0 0 0. 000111 000111
(iii) Modify the MATLAB function SMW3.m which uses SMW to solve the linear systems (An + UnVnT )x = bn, where
An =In, UnT = 1 1 ... 1 ,VT = 1 2 ... n ,b=[1,...,1]T andn=10k, 11...1 23...n+1
for k = 1 . . . 6. Measure the computation times.
Observation: How do the computation times for the solution of the linear systems scale with the dimension? Hint: make sure that you pass A as a sparse matrix to SMW.
If the inverses of A and A ̄ exist, then
see next page
Large Scale Optimization for Data Science, Jacek Gondzio 3 Problem 3
Define Your Personal Linear Program
Letmandddenotethemonthanddayofyourbirthday. Clearly,1≤m≤12and1≤d≤31. Defineanumber: α=100·m+d,andconsideranLP
min x1 + x2 s.t. x1 + 2x2 x2
x1
x1≥0 x2≥0
Please complete the routine
+ 6x3 + 2x3
x3 + 3x3 x3 ≥ 0
+ x4 + x4
− x4 x4≥0
+ 3x5 + 4x6 −x6
+x5+2x6 +3x5+ x6 x5≥0 x6≥0.
= zLP ≤3 ≤3 ≥α ≤ 2+α
function [ alpha, fval, x, reducedCost, activePrimalConsts] = solveLp(m, d)
to model this problem. Solve the problem using the provided MATLAB implementation of the primal- dual method for linear programming
function [x,y,s] = ipdriver(c,A,b,tol,maxit,printlevel)
The OUTPUT of the routine solveLp should produce answers to the following questions:
(i) What is your personal α (alpha)?
(ii) What is the optimal objective (fval)?
(iii) Whatarethevaluesofvariables xˆi,i=1,...,6(x)?
(iv) What are the corresponding dual slacks di, i = 1, . . . , 6 (reducedCost)?
(v) Which primal constraints are active (activePrimalConsts)?
Please return which indexes among i = 1, 2, 3, 4 correspond to active constraints. (Hint: this is a numerical method, it is unlikely that a variable will reach an exact zero value; instead, it is likely that it will be close to zero).
Problem 4
Study the available MATLAB implementation of the primal-dual method for linear programming.
function [x,y,s] = ipdriver(c,A,b,tol,maxit,printlevel)
Extend this algorithm in
function [x,y,s] = ipdriverQP(c,A,b,Q,tol,maxit,printlevel)
so that it solves quadratic programming problems.
Define Your Personal Quadratic Program
Consider an extension of your personal LP in which a quadratic term has been added to the objective:
giving the following QP: min zQP
z Q P = z L P + 1 x 21 + 5 x 25 22
+x4
s.t. x1 x1
x1 ≥ 0
Complete the routine
+ 2x2 x2
x2≥0
+ 2x3
x3 −
≤3 −x6 ≤3
≥ α ≤ 2+α
function [ alpha, fval, x] = solveQP(m, d)
x4 x3≥0 x4≥0
+
+3x5+ x6
+ 3x3
that uses your QP interior point algorithm to solve this problem.
The OUTPUT of your routine should produce answers to the following questions:
(i) What is your personal α (alpha)?
(ii) What is the optimal objective (fval)?
(iii) Whatarethevaluesofvariables xˆi,i=1,...,6(x)?
x5 + 2x6 x5≥0 x6≥0.