Linear regression

YALMIP started out as a parser for semidefinite programs, but is just as useful for linear and quadratic programming. As an example, we will use YALMIP to define a couple of regression problems.

Let us assume that we have data generated from a noisy linear regression y(t) = x(t)a+e(t). The goal is to estimate the parameter a, given y and x, and we will try 3 different approaches.

Create some data to work with.

a = [1 2 3 4 5 6]';
t = (0:0.02:2*pi)';
x = [sin(t) sin(2*t) sin(3*t) sin(4*t) sin(5*t) sin(6*t)];
y = x*a+(-4+8*rand(length(x),1));

Define the variable we are looking for

 a_hat = sdpvar(6,1);

By using a_hat and the regressors x, we can define the residuals (which also will be an sdpvar object, parameterized in a_hat)

residuals = y-x*a_hat;

To solve the L₁ regression problem (minimize sum of absolute values of residuals), we define a variable that will serve as a bound on the elements in |y-x*a_hat|.

bound = sdpvar(length(residuals),1);

Express that bound is larger than the absolute values of the residuals (note the simple definition of a double-sided constraint).

F = set(-bound < residuals < bound);

Minimize the sum of the bounds subject to the constraints in F.

solvesdp(F,sum(bound));

The optimal value is extracted using the overloaded double command.

a_L1 = double(a_hat);

The L₂ problem is easily solved as a QP problem without any constraints.

solvesdp([],residuals'*residuals);
a_L2 = double(a_hat);

YALMIP automatically detects that the objective is a convex quadratic function, and solves the problem using any installed QP solver. If no QP solver is found, the problem is converted to an SOCP, and if no dedicated SOCP solver exist, the SOCP is converted to an SDP.

Finally, we minimize the L_∞ norm. This corresponds to minimizing the largest (absolute value) residual. Introduce a scalar to bound the largest value in the vector residual (YALMIP uses MATLAB standard to compare scalars, vectors and matrices)

bound = sdpvar(1,1);
F = set(-bound < residuals < bound);

and minimize the bound.

solvesdp(F,bound);
a_Linf = double(a_hat);

Make sure to check out the nonlinear operator examples to see how you can simplify the code even more.