0

在 Haskell 中，岭回归可以表示为：

import Numeric.LinearAlgebra 

createReadout :: Matrix Double → Matrix Double → Matrix Double
createReadout a b = oA <\> oB
  where
   μ = 1e-4

   oA = (a <> (tr a)) + (μ * (ident $ rows a))
   oB = a <> (tr b)

但是，此操作非常耗费内存。这是一个简约的示例，它需要我的机器上超过 2GB 的空间，并且需要 3 分钟才能执行。

import Numeric.LinearAlgebra
import System.Random

createReadout :: Matrix Double -> Matrix Double -> Matrix Double
createReadout a b = oA <\> oB
  where
    mu = 1e-4
    oA = (a <> (tr a)) + (mu * (ident $ rows a))
    oB = a <> (tr b)

teacher :: [Int] -> Int -> Int -> Matrix Double
teacher labelsList cols' correctRow = fromBlocks $ f <$> labelsList
  where ones = konst 1.0 (1, cols')
        zeros = konst 0.0 (1, cols')
        rows' = length labelsList
        f i | i == correctRow = [ones]
            | otherwise = [zeros]

glue :: Element t => [Matrix t] -> Matrix t
glue xs = fromBlocks [xs]

main :: IO ()
main = do

  let n = 1500  -- <- The constant to be increased
      m = 10000
      cols' = 12

  g <- newStdGen

  -- Stub data
  let labels = take m . map (`mod` 10) . randoms $ g :: [Int]
      a = (n >< (cols' * m)) $ take (cols' * m * n) $ randoms g :: Matrix Double
      teachers = zipWith (teacher [0..9]) (repeat cols') labels
      b = glue teachers

  print $ maxElement $ createReadout a b
  return ()

$ cabal exec ghc -- -O2 Test.hs

$ time ./Test
./Test 190.16s user 5.22s system 106% cpu 3:03.93 total

问题是增加常数n，至少到 n = 4000，而 RAM 受到 5GB 的限制。理论上矩阵求逆运算需要的最小空间是多少？这种操作如何在空间方面进行优化？岭回归可以用更便宜的方法有效地替代吗？

4

2 回答 2

1

简单的 Gauss-Jordan 消元只占用存储输入和输出矩阵的空间加上常数辅助空间。如果我没看错，oA你需要反转的矩阵是nx n，所以这不是问题。

您的内存使用完全由存储输入矩阵支配，该矩阵a至少使用 1500 * 120000 * 8 = 1.34 GB。n = 4000将是 4000 * 120000 * 8 = 3.58 GB，这是您空间预算的一半以上。我不知道您正在使用什么矩阵库或它如何存储它的矩阵，但如果它们在 Haskell 堆上，那么 GC 效果很容易占空间使用量的另一个因素 2。

于 2016-12-03T15:51:24.187 回答
1

好吧，您可以摆脱 3*m + nxn 空间，但是我不确定这在数值上有多稳定。

基础是身份

inv( inv(Q) + A'*A)) = Q - Q*A'*R*A*Q
where R = inv( I + A*Q*A')

如果 A 是您的 A 矩阵并且

Q = inv( mu*I*mu*I) = I/(mu*mu)

那么岭回归的解决方案是

inv( inv(Q) + A'*A)) * A'*b

更多的代数展示

inv( inv(Q) + A'*A)) = (I - A'*inv( (mu2 + A*A'))*A)/mu2
where mu2 = mu*m

请注意，因为 A 是 nxm，所以 A*A' 是 nx n。

所以一种算法是

计算 C = A*A' + mu2

对C进行cholesky分解，即找到上三角U，使得U'*U = C

计算向量 y = A'*b

计算向量 z = A*y

求解 U'*u = z for u in z

求解 U*v = z for v in z

计算 w = A'*z

计算 x = (y - w)/mu2。

于 2016-12-06T18:04:43.967 回答