1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
|
{- |
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)])
--------------------------------------------------------------------------------
|
