import%20marimo%0A%0A__generated_with%20%3D%20%220.15.2%22%0Aapp%20%3D%20marimo.App()%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20!%5BMOSEK%20ApS%5D(https%3A%2F%2Fwww.mosek.com%2Fstatic%2Fimages%2Fbranding%2Fwebgraphmoseklogocolor.png%20)%0A%0A%20%20%20%20%23%20Linear%20Regression%20techniques%20using%20MOSEK%20Fusion%20API%0A%0A%0A%20%20%20%20Regression%20is%20one%20of%20the%20most%20used%20techniques%20in%20finance%2C%20statistic%2C%20biology%20and%20in%20general%20any%20application%20where%20it%20is%20necessary%20to%20construct%20models%20that%20approximate%20input%20data.%0A%0A%20%20%20%20The%20aim%20of%20this%20tutorial%20is%20two-fold%3A%0A%0A%20%20%20%201.%20to%20show%20some%20modeling%20techniques%20to%20define%20regression%20problems%3B%0A%20%20%20%202.%20to%20show%20how%20to%20efficiently%20implement%20those%20problems%20using%20MOSEK%20Fusion%20API.%0A%0A%0A%20%20%20%20This%20tutorial%20is%20largerly%20based%20on%3A%0A%0A%20%20%20%20%20%20%20%20%5B1%5D%20Schmelzer%2C%20T.%2C%20Hauser%2C%20R.%2C%20Andersen%2C%20E.%2C%20%26%20Dahl%2C%20J.%20(2013).%20Regression%20techniques%20for%20Portfolio%20Optimisation%20using%20MOSEK.%20arXiv%20preprint%20arXiv%3A1310.3397.%0A%20%20%20%20%20%20%20%20%5B2%5D%20Boyd%2C%20S.%2C%20%26%20Vandenberghe%2C%20L.%20(2004).%20Convex%20optimization.%20Cambridge%20university%20press.%0A%0A%0A%20%20%20%20The%20code%20is%20written%20in%20Python3%20and%20uses%20the%20%5BMOSEK%20Fusion%20API%5D(https%3A%2F%2Fmosek.com%2Fdocumentation).%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20mosek%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%2C%20SolutionStatus%0A%20%20%20%20import%20mosek.fusion.pythonic%0A%20%20%20%20print(%22MOSEK%20version%22%2C%20Model.getVersion())%0A%20%20%20%20return%20Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20Var%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%20The%20basics%20of%20linear%20regression%0A%0A%20%20%20%20A%20linear%20model%20is%20used%20to%20explain%20the%20relation%20between%20%24m%24%20input%20variables%20%24x_1%2C%5Cldots%2Cx_m%24%20and%20a%20dependent%20variable%20%24y%24%20(the%20output).%20To%20this%20extent%20we%20assume%20that%20there%20exists%20a%20set%20of%20weights%20%24w%5Cin%20%5Cmathbb%7BR%7D%5E%7Bm%7D%24%20such%20that%20ideally%0A%0A%20%20%20%20%24%24%0A%20%20%20%20y%20%3D%20x%5ET%20w%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24x%3D(x_1%2C%5Cldots%2Cx_m)%24.%0A%0A%20%20%20%20In%20practice%2C%20we%20are%20given%20%24n%24%20observations%20of%20the%20relation%20that%20links%20%24x%24%20and%20%24y%24.%20We%20store%20them%20row-wise%20in%20a%20matrix%20%24X%5Cin%20%5Cmathbb%7BR%7D%5E%7Bn%5Ctimes%20m%7D%24.%20Corresponding%20outputs%20%24y%3D(y_i%2C%5Cldots%2Cy_n)%24%20are%20known%20as%20well.%20If%20the%20relation%20we%20want%20to%20describe%20is%20truly%20linear%20we%20can%20hope%20to%20solve%20the%20linear%20system%0A%0A%20%20%20%20%24%24%0A%20%20%20%20X%20w%20%3Dy.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20But%20usually%20this%20is%20not%20the%20case%20and%20we%20can%20only%20settle%20for%20an%20approximation%20which%20minimizes%20(in%20some%20sense)%20the%20residual%20vector%0A%0A%20%20%20%20%24%24%0A%20%20%20%20r%20%3D%20Xw%20-%20y%2C%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20we%20obtain%20an%20optimization%20problem%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%5Cphi(r)%5C%5C%20%0A%20%20%20%20%20%20s.t.%20%26%20r%20%3D%20Xw%20-%20y%2C%20%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%26%20w%5Cin%20W.%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20choice%20of%20the%20cost%20function%20%24%5Cphi(%5Ccdot)%24%20determines%20the%20way%20errors%20on%20the%20residuals%20are%20weighted.%20Here%20we%20demonstrate%20the%20following%20variants%3A%0A%0A%20%20%20%20*%20%24%5Cell_2%24-norm%2C%0A%20%20%20%20*%20%24%5Cell_1%24-norm%2C%0A%20%20%20%20*%20%24%5Cell_p-%24norm%20%2C%0A%20%20%20%20*%20*Chebyshev*%20norm%20(also%20known%20as%20*Min-max*)%2C%0A%20%20%20%20*%20*Deadzone-linear*%2C%0A%20%20%20%20*%20*%24k%24-largest%20residuals*.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20create%20simple%20input%20data%20sets%20which%20attempt%20to%20lead%20to%20the%20solution%20%24w%3D(1%2C%5Cldots%2C1)%24%20and%20we%20add%20small%20normally%20distributed%20error%20as%20follows.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20sys%0A%20%20%20%20%23%20'%25matplotlib%20inline'%20command%20supported%20automatically%20in%20marimo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%0A%20%20%20%20def%20dataset(m%2C%20n)%3A%0A%20%20%20%20%20%20%20%20X%20%3D%20np.random.uniform(-1.0%2C%201.0%2C%20(n%2Cm))%0A%20%20%20%20%20%20%20%20y%20%3D%20X.dot(np.ones(m))%20%2B%20np.random.normal(0.0%2C%200.2%2C%20n)%0A%20%20%20%20%20%20%20%20return%20X%2C%20y%0A%20%20%20%20return%20dataset%2C%20np%2C%20plt%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20%24%5Cell_2%24-norm%0A%0A%20%20%20%20The%20regression%20problem%20with%20respect%20to%20the%20%24%5Cell_2%24-norm%20(i.e.%20linear%20least%20squares)%20has%20the%20conic%20quadratic%20formulation%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20t%5C%5C%0A%20%20%20%20%5Ctext%7Bs.t.%7D%20%26%20t%5Cgeq%20%5C%7CXw-y%5C%7C_2%20%5Cquad%20%5Cmathrm%7Bi.e.%7D%5C%20(t%2CXw-y)%5Cin%20%5Cmathcal%7BQ%7D%5E%7Bn%2B1%7D%2C%5C%5C%0A%20%20%20%20%20%26%20t%5Cin%5Cmathbb%7BR%7D%2C%20w%5Cin%5Cmathbb%7BR%7D%5Em.%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20This%20model%20models%20the%20maximum%20likelihood%20of%20the%20fit%20precisely%20under%20the%20assumption%20that%20residuals%20have%20normal%20distribution.%0A%0A%20%20%20%20Below%20is%20the%20model%20in%20Fusion.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20dataset)%3A%0A%20%20%20%20def%20l2norm(X%2Cy)%3A%0A%20%20%20%20%20%20%20%20n%2Cm%3DX.shape%0A%0A%20%20%20%20%20%20%20%20M%20%3D%20Model(%22l2norm%22)%0A%20%20%20%20%20%20%20%20w%20%3D%20M.variable(m)%0A%20%20%20%20%20%20%20%20t%20%3D%20M.variable()%0A%0A%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20t)%0A%0A%20%20%20%20%20%20%20%20res%20%3D%20X%20%40%20w%20-%20y%0A%20%20%20%20%20%20%20%20M.constraint(Expr.vstack(t%2C%20res)%2C%20Domain.inQCone())%0A%0A%20%20%20%20%20%20%20%20%23M.setLogHandler(sys.stdout)%0A%20%20%20%20%20%20%20%20M.solve()%0A%0A%20%20%20%20%20%20%20%20%23%20Return%20the%20weights%20vector%20and%20the%20residuals%0A%20%20%20%20%20%20%20%20w%20%3D%20w.level()%0A%20%20%20%20%20%20%20%20return%20w%2C%20X.dot(w)-y%0A%0A%20%20%20%20X%2Cy%20%3D%20dataset(10%2C1000)%0A%20%20%20%20w%2C%20res%20%3D%20l2norm(X%2Cy)%0A%20%20%20%20return%20X%2C%20res%2C%20w%2C%20y%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22If%20everything%20went%20well%2C%20we%20should%20have%20recovered%20the%20weight%20vector%20of%20all%20ones%20and%20normally%20distributed%20residuals%2C%20as%20below%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np%2C%20plt%2C%20res%2C%20w)%3A%0A%20%20%20%20print(%22%7Cw-1%7C_2%3A%22%2C%20np.linalg.norm(w-np.ones(10)))%0A%20%20%20%20plt.figure(figsize%3D(6%2C4))%20%0A%20%20%20%20plt.hist(res%2C40)%0A%20%20%20%20plt.title(%22Histogram%20of%20residuals%20(l2%20norm)%22)%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20%24%5Cell_1%24-norm%0A%0A%20%20%20%20Under%20the%20%24%5Cell_1%24-norm%20we%20minimize%20%24%5Csum_i%20%7Cr_i%7C%24.%20The%20corresponding%20linear%20formulation%20is%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%20%26%20%5Csum_i%20t_i%20%5C%5C%0A%20%20%20%20%5Ctext%7Bs.t.%7D%20%26%20t%5Cgeq%20Xw-y%2C%5C%5C%0A%20%20%20%20%20%26%20t%20%5Cgeq%20-(Xw-y)%2C%20%5C%5C%0A%20%20%20%20%20%26%20w%5Cin%5Cmathbb%7BR%7D%5Em%2C%20t%5Cin%5Cmathbb%7BR%7D%5En%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%20%0A%0A%20%20%20%20which%20simply%20expresses%20the%20bound%20%24t_i%5Cgeq%20%7Cr_i%7C%24.%20Below%20is%20a%20corresponding%20Fusion%20model.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Expr%2C%20Model%2C%20ObjectiveSense%2C%20X%2C%20np%2C%20plt%2C%20y)%3A%0A%20%20%20%20def%20l1norm(X%2C%20y)%3A%0A%20%20%20%20%20%20%20%20n%2C%20m%20%3D%20X.shape%0A%20%20%20%20%20%20%20%20M%20%3D%20Model('l1norm')%0A%20%20%20%20%20%20%20%20w%20%3D%20M.variable(m)%0A%20%20%20%20%20%20%20%20t%20%3D%20M.variable(n)%0A%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20Expr.sum(t))%0A%20%20%20%20%20%20%20%20res%20%3D%20X%20%40%20w%20-%20y%0A%20%20%20%20%20%20%20%20M.constraint(t%20%3E%3D%20res)%0A%20%20%20%20%20%20%20%20M.constraint(t%20%3E%3D%20-res)%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20w%20%3D%20w.level()%0A%20%20%20%20%20%20%20%20return%20(w%2C%20X.dot(w)%20-%20y)%0A%20%20%20%20w_1%2C%20res_1%20%3D%20l1norm(X%2C%20y)%0A%20%20%20%20print('%7Cw-1%7C_2%3A'%2C%20np.linalg.norm(w_1%20-%20np.ones(10)))%0A%20%20%20%20plt.hist(res_1%2C%2040)%0A%20%20%20%20plt.title('Histogram%20of%20residuals%20(l1%20norm)')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20%24%5Cell_p%24-norm%0A%0A%20%20%20%20The%20%24%5Cell_p%24%20norm%0A%0A%20%20%20%20%24%24%5Cphi(r)%20%3D%20%5C%7C%20r%20%5C%7C_p%20%3D%20%5Cleft(%20%5Csum_i%20%7Cr_i%7C%5Ep%20%5Cright)%5E%7B1%2Fp%7D%24%24%0A%0A%20%20%20%20For%20a%20general%20%24p%3E1%24%20we%20can%20model%20it%20using%20the%20power%20cone.%20The%20standard%20representation%20of%20the%20epigraph%20%24t%5Cgeq%5C%7Cr%5C%7C_p%24%20as%20a%20conic%20model%20can%20be%20found%20in%20https%3A%2F%2Fdocs.mosek.com%2Fmodeling-cookbook%2Fpowo.html%23norm-cones%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bl%7D%0A%20%20%20%20%5Csum_i%20s_i%20%3D%20t%5C%5C%0A%20%20%20%20s_it%5E%7Bp-1%7D%20%5Cgeq%20%7Cr_i%7C%5Ep%2C%20%5C%20i%3D1%2C%5Cldots%2Cn%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20and%20the%20last%20constraint%20can%20be%20expressed%20using%20a%20power%20cone.%20In%20Fusion%20all%20of%20these%20individual%20conic%20constraints%20can%20be%20specified%20at%20once%20in%20matrix%20for%2C%20where%2C%20as%20always%20in%20Fusion%2C%20%22matrix%20in%20a%20cone%22%20is%20a%20shorthand%20for%20%22each%20row%20of%20the%20matrix%20is%20in%20the%20cone%22.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20Var%2C%20X%2C%20np%2C%20plt%2C%20y)%3A%0A%20%20%20%20def%20lpnorm(X%2C%20y%2C%20p)%3A%0A%20%20%20%20%20%20%20%20n%2C%20m%20%3D%20X.shape%0A%20%20%20%20%20%20%20%20M%20%3D%20Model('lpnorm')%0A%20%20%20%20%20%20%20%20w%20%3D%20M.variable(m)%0A%20%20%20%20%20%20%20%20t%20%3D%20M.variable()%0A%20%20%20%20%20%20%20%20s%20%3D%20M.variable(n)%0A%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20t)%0A%20%20%20%20%20%20%20%20M.constraint(t%20%3D%3D%20Expr.sum(s))%0A%20%20%20%20%20%20%20%20res%20%3D%20X%20%40%20w%20-%20y%0A%20%20%20%20%20%20%20%20M.constraint(Expr.hstack(s%2C%20Var.repeat(t%2C%20n)%2C%20res)%2C%20Domain.inPPowerCone(1.0%20%2F%20p))%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20w%20%3D%20w.level()%0A%20%20%20%20%20%20%20%20return%20(w%2C%20X.dot(w)%20-%20y)%0A%20%20%20%20w_2%2C%20res_2%20%3D%20lpnorm(X%2C%20y%2C%20p%3D2)%0A%20%20%20%20print('%7Cw-1%7C_2%3A'%2C%20np.linalg.norm(w_2%20-%20np.ones(10)))%0A%20%20%20%20plt.hist(res_2%2C%2040)%0A%20%20%20%20plt.title('Histogram%20of%20residuals%20(lp%20norm)')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%20Deadzone-linear%0A%0A%20%20%20%20This%20approach%20steams%20from%20the%20%24%5Cell_1%24%20norm%20approximation%2C%20but%20ignoring%20small%20errors.%20Given%20a%20threshold%20%24a%5Cgeq%200%24%2C%20we%20define%20%0A%0A%20%20%20%20%24%24%20%5Cphi(r)%20%3D%20%5Csum_i%20%5Cmax(0%2C%20%7Cr_i%7C-a%20)%24%24%0A%0A%20%20%20%20The%20conic%20optimization%20model%20takes%20the%20form%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%20%5Csum_i%20t_i%5C%5C%0A%20%20%20%20%5Ctext%7Bs.t.%7D%20%26%20t%20%5Cgeq%20%7CXw-y%7C-a%2C%5C%5C%0A%20%20%20%20%26%20t%5Cgeq%200%20%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%20%0A%0A%20%20%20%20Again%2C%20removing%20the%20absolute%20value%20we%20obtain%20in%20compact%20form%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20%20%5Csum_i%20t_i%5C%5C%0A%20%20%20%20%5Ctext%7Bs.t.%7D%20%26%20t%20%5Cgeq%20Xw-y-a%2C%5C%5C%0A%20%20%20%20%20%20%26%20t%20%5Cgeq%20-(Xw-y)-a%2C%5C%5C%0A%20%20%20%20%26%20t%5Cgeq%200.%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20X%2C%20np%2C%20plt%2C%20y)%3A%0A%20%20%20%20def%20deadzone(X%2C%20y%2C%20a)%3A%0A%20%20%20%20%20%20%20%20n%2C%20m%20%3D%20X.shape%0A%20%20%20%20%20%20%20%20M%20%3D%20Model('deadzone')%0A%20%20%20%20%20%20%20%20w%20%3D%20M.variable(m)%0A%20%20%20%20%20%20%20%20t%20%3D%20M.variable(n%2C%20Domain.greaterThan(0))%0A%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20Expr.sum(t))%0A%20%20%20%20%20%20%20%20res%20%3D%20X%20%40%20w%20-%20y%0A%20%20%20%20%20%20%20%20M.constraint(t%20%3E%3D%20res%20-%20a)%0A%20%20%20%20%20%20%20%20M.constraint(t%20%3E%3D%20-res%20-%20a)%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20w%20%3D%20w.level()%0A%20%20%20%20%20%20%20%20return%20(w%2C%20X.dot(w)%20-%20y)%0A%20%20%20%20w_3%2C%20res_3%20%3D%20deadzone(X%2C%20y%2C%200.05)%0A%20%20%20%20print('%7Cw-1%7C_2%3A'%2C%20np.linalg.norm(w_3%20-%20np.ones(10)))%0A%20%20%20%20plt.hist(res_3%2C%2040)%0A%20%20%20%20plt.title('Histogram%20of%20residuals%20(deadzone)')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20Chebyshev%20Approximation%0A%0A%20%20%20%20Sometimes%20referred%20as%20the%20*min-max*%20approximation%2C%20it%20uses%20%0A%0A%20%20%20%20%24%24%5Cphi(r)%20%3D%20%5Cmax_i(%20%7Cr_i%7C%20)%20%3D%20%5C%7Cr%5C%7C_%7B%5Cinfty%7D%24%24%0A%0A%20%20%20%20as%20penalty%20function.%20It%20is%20a%20limit%20case%20of%20%24%5Cell_p-%24norms%20as%20%24p%5Crightarrow%20%5Cinfty%24.%20The%20effect%20is%20to%20only%20look%20at%20the%20largest%20error%20in%20absolute%20value.%20The%20formulation%20of%20our%20problem%20is%20therefore%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20t%5C%5C%0A%20%20%20%20%5Ctext%7Bs.t.%7D%20%26%20te%20%5Cgeq%20Xw-y%2C%5C%5C%0A%20%20%20%20%20%20%26%20te%20%5Cgeq%20-(Xw-y)%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24e%24%20is%20the%20all-ones%20vector.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Model%2C%20ObjectiveSense%2C%20Var%2C%20X%2C%20np%2C%20plt%2C%20y)%3A%0A%20%20%20%20def%20minmax(X%2C%20y)%3A%0A%20%20%20%20%20%20%20%20n%2C%20m%20%3D%20X.shape%0A%20%20%20%20%20%20%20%20M%20%3D%20Model('minmax')%0A%20%20%20%20%20%20%20%20w%20%3D%20M.variable(m)%0A%20%20%20%20%20%20%20%20t%20%3D%20M.variable()%0A%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20t)%0A%20%20%20%20%20%20%20%20res%20%3D%20X%20%40%20w%20-%20y%0A%20%20%20%20%20%20%20%20t_repeat%20%3D%20Var.repeat(t%2C%20n)%0A%20%20%20%20%20%20%20%20M.constraint(t_repeat%20%3E%3D%20res)%0A%20%20%20%20%20%20%20%20M.constraint(t_repeat%20%3E%3D%20-res)%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20w%20%3D%20w.level()%0A%20%20%20%20%20%20%20%20return%20(w%2C%20X.dot(w)%20-%20y)%0A%20%20%20%20w_4%2C%20res_4%20%3D%20minmax(X%2C%20y)%0A%20%20%20%20print('%7Cw-1%7C_2%3A'%2C%20np.linalg.norm(w_4%20-%20np.ones(10)))%0A%20%20%20%20plt.hist(res_4%2C%2040)%0A%20%20%20%20plt.title('Histogram%20of%20residuals%20(minmax)')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20Sum%20of%20the%20largest%20%24k%24%20residuals%0A%0A%20%20%20%20The%20penalty%20function%20%24%5Cphi%24%20is%20defined%20as%20%0A%0A%20%20%20%20%24%24%5Cphi(u)%20%3D%20%5Csum_%7Bi%3D1%7D%5Ek%20%7C%20r_%7B%5Bi%5D%7D%20%7C%2C%20%24%24%0A%0A%20%20%20%20where%20the%20notation%20%24r_%7B%5Bi%5D%7D%24%20indicates%20the%20%24i%24-th%20element%20of%20%24r%24%20sorted.%20It%20means%20that%20the%20penalty%20depends%20on%20the%20%24k%24%20largest%20(in%20absolute%20value)%20residuals.%20This%20is%20a%20generalization%20of%20the%20minmax%20norm%2C%20which%20here%20corresponds%20to%20%24k%3D1%24%2C%20and%20of%20the%20%24%5Cell_1%24-norm%2C%20which%20corresponds%20to%20%24k%3Dn%24.%20As%20shown%20in%20%5Bin%20our%20documentation%20of%20sum%20of%20largest%20elements%5D(https%3A%2F%2Fdocs.mosek.com%2Fmodeling-cookbook%2Flinear.html%23sum-of-largest-elements)%20this%20problem%20can%20be%20modeled%20in%20a%20linear%20fashion%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20ks%2B%5Csum%20u_i%5C%5C%0A%20%20%20%20%5Ctext%7Bs.t.%7D%20%26%20se%20%2B%20u%20%5Cgeq%20Xw%20-y%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%26%20se%20%2B%20u%20%5Cgeq%20-(Xw%20-y)%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%26%20u%5Cgeq%200%2C%5C%5C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%26%20u%5Cin%20%5Cmathbb%7BR%7D%5En%2Cs%5Cin%20%5Cmathbb%7BR%7D%2Cw%5Cin%5Cmathbb%7BR%7D%5Em.%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20Var%2C%20X%2C%20np%2C%20plt%2C%20y)%3A%0A%20%20%20%20def%20sumkmax(X%2C%20y%2C%20k)%3A%0A%20%20%20%20%20%20%20%20n%2C%20m%20%3D%20X.shape%0A%20%20%20%20%20%20%20%20M%20%3D%20Model('sumkmax')%0A%20%20%20%20%20%20%20%20w%20%3D%20M.variable(m)%0A%20%20%20%20%20%20%20%20s%20%3D%20M.variable()%0A%20%20%20%20%20%20%20%20u%20%3D%20M.variable(n%2C%20Domain.greaterThan(0))%0A%20%20%20%20%20%20%20%20M.objective(ObjectiveSense.Minimize%2C%20Expr.sum(u)%20%2B%20k%20*%20s)%0A%20%20%20%20%20%20%20%20res%20%3D%20X%20%40%20w%20-%20y%0A%20%20%20%20%20%20%20%20lhs%20%3D%20u%20%2B%20Var.repeat(s%2C%20n)%0A%20%20%20%20%20%20%20%20M.constraint(lhs%20%3E%3D%20res)%0A%20%20%20%20%20%20%20%20M.constraint(lhs%20%3E%3D%20-res)%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20w%20%3D%20w.level()%0A%20%20%20%20%20%20%20%20return%20(w%2C%20X.dot(w)%20-%20y)%0A%20%20%20%20w_5%2C%20res_5%20%3D%20sumkmax(X%2C%20y%2C%2020)%0A%20%20%20%20print('%7Cw-1%7C_2%3A'%2C%20np.linalg.norm(w_5%20-%20np.ones(10)))%0A%20%20%20%20plt.hist(res_5%2C%2040)%0A%20%20%20%20plt.title('Histogram%20of%20residuals%20(sumkmax)')%0A%20%20%20%20plt.show()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%23%20Some%20observations%0A%0A%20%20%20%20All%20the%20regression%20models%20presented%20in%20this%20tutorial%20lead%20to%20either%20a%20linear%20or%20a%20second-order%20cone%20or%20a%20power%20cone%20%20optimization%20problem.%20That%20means%20that%0A%0A%20%20%20%20*%20they%20can%20be%20solved%20in%20polynomial%20time%20by%20standard%20optimization%20algorithms%3B%0A%20%20%20%20*%20add%20additional%20(conic)constraints%20on%20the%20weights%20can%20be%20easily%20included%3B%0A%20%20%20%20*%20sparsity%20in%20the%20input%2C%20i.e.%20in%20%24X%24%2C%20can%20be%20easily%20handled%20and%20usually%20leads%20to%20great%20savings%20in%20time%20and%20memory%3B%0A%20%20%20%20*%20the%20solver%20will%20automatically%20choose%20between%20the%20primal%20and%20dual%20formulation%3B%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3E%3Cimg%20alt%3D%22Creative%20Commons%20License%22%20style%3D%22border-width%3A0%22%20src%3D%22https%3A%2F%2Fi.creativecommons.org%2Fl%2Fby%2F4.0%2F80x15.png%22%20%2F%3E%3C%2Fa%3E%3Cbr%20%2F%3EThis%20work%20is%20licensed%20under%20a%20%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3ECreative%20Commons%20Attribution%204.0%20International%20License%3C%2Fa%3E.%20The%20**MOSEK**%20logo%20and%20name%20are%20trademarks%20of%20%3Ca%20href%3D%22http%3A%2F%2Fmosek.com%22%3EMosek%20ApS%3C%2Fa%3E.%20The%20code%20is%20provided%20as-is.%20Compatibility%20with%20future%20release%20of%20**MOSEK**%20or%20the%20%60Fusion%20API%60%20are%20not%20guaranteed.%20For%20more%20information%20contact%20our%20%5Bsupport%5D(mailto%3Asupport%40mosek.com).%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20return%20(mo%2C)%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
c3388793d262d93161525f30fdd656aef23a54d2c12a1d741676a05467d5398b