{- | Module : Numeric.LinearProgramming.L1 Copyright : (c) Alberto Ruiz 2011-14 Stability : provisional Linear system solvers in the L_1 norm using linear programming. -} ----------------------------------------------------------------------------- module Numeric.LinearProgramming.L1 ( l1Solve, l1SolveGT, l1SolveO, lInfSolveO, l1SolveU, ) where import Numeric.LinearAlgebra import Numeric.LinearProgramming -- | L_inf solution of overconstrained system Ax=b. -- -- @argmin_x ||Ax-b||_inf@ lInfSolveO :: Matrix Double -> Vector Double -> Vector Double lInfSolveO a b = fromList (take n x) where n = cols a as = toRows a bs = toList b c1 = zipWith (mk (1)) as bs c2 = zipWith (mk (-1)) as bs mk sign a_i b_i = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(n+1)]) :<=: (sign * b_i) p = Sparse (c1++c2) Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ [1])) p (map Free [1..(n+1)]) -------------------------------------------------------------------------------- -- | L_1 solution of overconstrained system Ax=b. -- -- @argmin_x ||Ax-b||_1@ l1SolveO :: Matrix Double -> Vector Double -> Vector Double l1SolveO a b = fromList (take n x) where n = cols a m = rows a as = toRows a bs = toList b ks = [1..] c1 = zipWith3 (mk (1)) as bs ks c2 = zipWith3 (mk (-1)) as bs ks mk sign a_i b_i k = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(k+n)]) :<=: (sign * b_i) p = Sparse (c1++c2) Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ replicate m 1)) p (map Free [1..(n+m)]) -------------------------------------------------------------------------------- -- | L1 solution of underconstrained linear system Ax=b. -- -- @argmin_x ||x||_1 such that Ax=b@ l1SolveU :: Matrix Double -> Vector Double -> Vector Double l1SolveU a y = fromList (take n x) where n = cols a c1 = map (\k -> [ 1#k, -1#k+n] :<=: 0) [1..n] c2 = map (\k -> [-1#k, -1#k+n] :<=: 0) [1..n] c3 = zipWith (:==:) (map sp \$ toRows a) (toList y) sp v = zipWith (#) (toList v) [1..] p = Sparse (c1 ++ c2 ++ c3) Optimal (_j,x) = simplex (Minimize (replicate n 0 ++ replicate n 1)) p (map Free [1..(2*n)]) -------------------------------------------------------------------------------- -- | Solution in the L_1 norm, with L_1 regularization, of a linear system @Ax=b@. -- -- @argmin_x λ||x||_1 + ||Ax-b||_1@ l1Solve :: Double -- ^ λ -> Matrix Double -- ^ A -> Vector Double -- ^ b -> Vector Double -- ^ x l1Solve λ a b = fromList (take n x) where n = cols a m = rows a as = toRows a bs = toList b c1Res = zipWith3 (mkR (1)) as bs [1..m] c2Res = zipWith3 (mkR (-1)) as bs [1..m] mkR sign a_i b_i k = (zipWith (#) (toList (scale sign a_i)) [1..] ++ [-1#(k+2*n)]) :<=: (sign * b_i) c1Sol = map (\k -> [ 1#k, -1#k+n] :<=: 0) [1..n] c2Sol = map (\k -> [-1#k, -1#k+n] :<=: 0) [1..n] p = Sparse (c1Res++c2Res++c1Sol++c2Sol) cost = replicate n 0 ++ replicate n λ ++ replicate m 1 Optimal (_j,x) = simplex (Minimize cost) p (map Free [1..(2*n+m)]) -------------------------------------------------------------------------------- -- | Solution in the L_1 norm, with L_1 regularization, of a system of linear inequalities @Ax>=b@. -- -- @argmin_x λ||x||_1 + ||step(b-Ax)||_1@ l1SolveGT :: Double -- ^ λ -> Matrix Double -- ^ A -> Vector Double -- ^ b -> Vector Double -- ^ x l1SolveGT λ a b = fromList (take n x) where n = cols a m = rows a as = toRows a bs = toList b cRes = zipWith3 mkR as bs [1..m] mkR a_i b_i k = (zipWith (#) (toList a_i) [1..] ++ [1#(k+2*n)]) :>=: (b_i) c1Sol = map (\k -> [ 1#k, -1#k+n] :<=: 0) [1..n] c2Sol = map (\k -> [-1#k, -1#k+n] :<=: 0) [1..n] p = Sparse (cRes++c1Sol++c2Sol) cost = replicate n 0 ++ replicate n λ ++ replicate m 1 Optimal (_j,x) = simplex (Minimize cost) p (map Free [1..(2*n)]) --------------------------------------------------------------------------------