岭回归需要多少space？

Question

在Haskell中，ridge regression可以表示为：

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，而内存限制在5GB。理论上矩阵求逆运算需要的最小space是多少？如何根据 space 优化此操作？可以用更便宜的方法有效地替代岭回归吗？

Answer 1

简单的Gauss-Jordan消元只需要space来存储输入输出矩阵加上常量辅助space。如果我没看错的话，你需要反转的矩阵 oA 是 n x n 所以这不是问题。

您的内存使用量完全由存储输入矩阵 a 支配，它至少使用 1500 * 120000 * 8 = 1.34 GB。 n = 4000 将是 4000 * 120000 * 8 = 3.58 GB，这超过了您 space 预算的一半。我不知道你使用的是什么矩阵库，也不知道它是如何存储矩阵的，但如果它们在 Haskell 堆上，那么 GC 效应可以很容易地解释 space 使用中的另一个因素 2。

Answer 2

好吧，你可以用 3*m + nxn space，但我不确定这在数值上有多稳定。

基础是身份

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 是 n x m，所以 A*A' 是 n x 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。

岭回归需要多少space？

How much space does ridge regression require?

haskell

machine-learning

linear-algebra

hmatrix