blob: 2e8f94aa8c4465a04c9979582c07d1f4b2d2ffbf (
plain)
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
|
{- |
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 (
l1SolveO, lInfSolveO,
l1SolveU,
) where
import Numeric.LinearAlgebra
import Numeric.LinearProgramming
-- | L_Inf solution of overconstrained system Ax=b.
--
-- Find 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.
--
-- Find 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.
--
-- Find 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)])
|
