[37] | 1 | <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN"> |
---|
| 2 | <html> |
---|
| 3 | |
---|
| 4 | <head> |
---|
| 5 | <meta http-equiv="Content-Language" content="en-us"> |
---|
| 6 | <title>YALMIP Example : Linear regression</title> |
---|
| 7 | <meta http-equiv="Content-Type" content="text/html; charset=windows-1251"> |
---|
| 8 | <meta content="Microsoft FrontPage 6.0" name="GENERATOR"> |
---|
| 9 | <meta name="ProgId" content="FrontPage.Editor.Document"> |
---|
| 10 | <link href="yalmip.css" type="text/css" rel="stylesheet"> |
---|
| 11 | <base target="_self"> |
---|
| 12 | </head> |
---|
| 13 | |
---|
| 14 | <body leftMargin="0" topMargin="0"> |
---|
| 15 | |
---|
| 16 | <div align="left"> |
---|
| 17 | <table border="0" cellpadding="4" cellspacing="3" style="border-collapse: collapse" bordercolor="#000000" width="100%" align="left" height="100%"> |
---|
| 18 | <tr> |
---|
| 19 | <td width="100%" align="left" height="100%" valign="top"> |
---|
| 20 | <h2>Linear regression</h2> |
---|
| 21 | <hr noShade SIZE="1" color="#000000"> |
---|
| 22 | <p>YALMIP started out as a parser for semidefinite programs, but is |
---|
| 23 | just as useful for linear and quadratic programming. As an example, |
---|
| 24 | we will use YALMIP to define a couple of regression problems.</p> |
---|
| 25 | <p>Let us assume that we have data generated from a noisy linear |
---|
| 26 | regression <strong>y(t) = x(t)a+e(t).</strong> The goal is to |
---|
| 27 | estimate the parameter <strong>a</strong>, given <strong>y</strong> |
---|
| 28 | and <strong>x</strong>, and we will try 3 different approaches.</p> |
---|
| 29 | <p>Create some data to work with.</p> |
---|
| 30 | <table cellPadding="10" width="100%"> |
---|
| 31 | <tr> |
---|
| 32 | <td class="xmpcode"> |
---|
| 33 | <pre>a = [1 2 3 4 5 6]'; |
---|
| 34 | t = (0:0.02:2*pi)'; |
---|
| 35 | x = [sin(t) sin(2*t) sin(3*t) sin(4*t) sin(5*t) sin(6*t)]; |
---|
| 36 | y = x*a+(-4+8*rand(length(x),1));</pre> |
---|
| 37 | </td> |
---|
| 38 | </tr> |
---|
| 39 | </table> |
---|
| 40 | <p>Define the variable we are looking for</p> |
---|
| 41 | <table cellPadding="10" width="100%"> |
---|
| 42 | <tr> |
---|
| 43 | <td class="xmpcode"> |
---|
| 44 | <pre> a_hat = sdpvar(6,1);</pre> |
---|
| 45 | </td> |
---|
| 46 | </tr> |
---|
| 47 | </table> |
---|
| 48 | <p>By using <code>a_hat</code> and the regressors |
---|
| 49 | <font face="Courier New" color="#0000c0">x</font>, we can define the |
---|
| 50 | residuals (which also will be an <a href="reference.htm#sdpvar">sdpvar</a> |
---|
| 51 | object, parameterized in <font face="Courier New" color="#0000c0"> |
---|
| 52 | a_hat</font>)</p> |
---|
| 53 | <table cellPadding="10" width="100%"> |
---|
| 54 | <tr> |
---|
| 55 | <td class="xmpcode"> |
---|
| 56 | <pre>residuals = y-x*a_hat;</pre> |
---|
| 57 | </td> |
---|
| 58 | </tr> |
---|
| 59 | </table> |
---|
| 60 | <p>To solve the L<sub>1</sub> regression problem (minimize sum of absolute |
---|
| 61 | values of residuals), we define a variable that will serve as a bound |
---|
| 62 | on the elements in |<font face="Courier New" color="#0000c0">y-x*a_hat</font>|.</p> |
---|
| 63 | <table cellPadding="10" width="100%"> |
---|
| 64 | <tr> |
---|
| 65 | <td class="xmpcode"> |
---|
| 66 | <pre>bound = sdpvar(length(residuals),1);</pre> |
---|
| 67 | </td> |
---|
| 68 | </tr> |
---|
| 69 | </table> |
---|
| 70 | <p>Express that <font color="#0000c0" face="Courier New">bound</font> |
---|
| 71 | is larger than the absolute values of the residuals (note the simple |
---|
| 72 | definition of a double-sided constraint).</p> |
---|
| 73 | <table cellPadding="10" width="100%"> |
---|
| 74 | <tr> |
---|
| 75 | <td class="xmpcode"> |
---|
| 76 | <pre>F = set(-bound < residuals < bound);</pre> |
---|
| 77 | </td> |
---|
| 78 | </tr> |
---|
| 79 | </table> |
---|
| 80 | <p>Minimize the sum of the bounds subject to the constraints in |
---|
| 81 | <font face="Courier New" color="#0000c0">F</font>.</p> |
---|
| 82 | <table cellPadding="10" width="100%"> |
---|
| 83 | <tr> |
---|
| 84 | <td class="xmpcode"> |
---|
| 85 | <pre>solvesdp(F,sum(bound));</pre> |
---|
| 86 | </td> |
---|
| 87 | </tr> |
---|
| 88 | </table> |
---|
| 89 | <p>The optimal value is extracted using the overloaded |
---|
| 90 | <a href="reference.htm#double">double</a> |
---|
| 91 | command.</p> |
---|
| 92 | <table cellPadding="10" width="100%"> |
---|
| 93 | <tr> |
---|
| 94 | <td class="xmpcode"> |
---|
| 95 | <pre>a_L1 = double(a_hat);</pre> |
---|
| 96 | </td> |
---|
| 97 | </tr> |
---|
| 98 | </table> |
---|
| 99 | <p>The L<sub>2</sub> problem is easily solved as a QP problem without any |
---|
| 100 | constraints.</p> |
---|
| 101 | <table cellPadding="10" width="100%"> |
---|
| 102 | <tr> |
---|
| 103 | <td class="xmpcode"> |
---|
| 104 | <pre>solvesdp([],residuals'*residuals); |
---|
| 105 | a_L2 = double(a_hat);</pre> |
---|
| 106 | </td> |
---|
| 107 | </tr> |
---|
| 108 | </table> |
---|
| 109 | <p>YALMIP automatically detects that the objective is a convex |
---|
| 110 | quadratic function, and solves the problem using any installed QP |
---|
| 111 | solver. If no QP solver is found, the problem is converted to an |
---|
| 112 | SOCP, and if no dedicated SOCP solver exist, the SOCP is converted to |
---|
| 113 | an SDP. </p> |
---|
| 114 | <p>Finally, we minimize the L<sub>∞</sub> norm. This corresponds to minimizing |
---|
| 115 | the largest (absolute value) residual. Introduce a scalar to bound |
---|
| 116 | the largest value in the vector <font color="#0000c0" face="Courier New"> |
---|
| 117 | residual</font> |
---|
| 118 | (YALMIP uses MATLAB standard to compare scalars, vectors and |
---|
| 119 | matrices)</p> |
---|
| 120 | <table cellPadding="10" width="100%"> |
---|
| 121 | <tr> |
---|
| 122 | <td class="xmpcode"> |
---|
| 123 | <pre>bound = sdpvar(1,1); |
---|
| 124 | F = set(-bound < residuals < bound);</pre> |
---|
| 125 | </td> |
---|
| 126 | </tr> |
---|
| 127 | </table> |
---|
| 128 | <p>and minimize the bound.</p> |
---|
| 129 | <table cellPadding="10" width="100%"> |
---|
| 130 | <tr> |
---|
| 131 | <td class="xmpcode"> |
---|
| 132 | <pre>solvesdp(F,bound); |
---|
| 133 | a_Linf = double(a_hat);</pre> |
---|
| 134 | </td> |
---|
| 135 | </tr> |
---|
| 136 | </table> |
---|
| 137 | <p><img border="0" src="demoicon.gif" width="16" height="16">Make sure to check out the <a href="extoperators.htm"> |
---|
| 138 | nonlinear operator examples</a> to see how you can |
---|
| 139 | simplify |
---|
| 140 | the code even more.</td> |
---|
| 141 | </tr> |
---|
| 142 | </table> |
---|
| 143 | </div> |
---|
| 144 | |
---|
| 145 | </body> |
---|
| 146 | |
---|
| 147 | </html> |
---|