Idiomatic Vector Algebra in Haskell

Question

As a way to practice with the vector library in Haskell, I'm trying to rewrite a Nelder-Mead minimization algorithm I had previously written in C. So far I've been having a bit of trouble translating some vector operations idiomatically.

For instance, consider a function which finds the centroid of n vectors out of a list of n+1 (filtering away one index),

In C, this can be written as

static void get_centroid(double **s, int n, int iz,
                         double *C)
{
  for (int i = 0; i < n+1; i++) {
    if (i != iz) {
      for (int j = 0; j < n; j++)
        C[j] += s[i][j];
    }
  }
  for (int j = 0; j < n; j++)
    C[j] /= n;
}

I tried translating this into Haskell, and ended up with the following

import Data.Vector
import qualified Data.Vector as V

type Node = Vector Double

type Simplex = Vector Node

centroid :: Simplex -> Int -> Node
centroid s iz = V.map (/ (fromIntegral $ V.length s)) $ V.zipWith (-) v (s ! iz) 
    where v = V.foldl go V.empty s 
        where go a b = V.zipWith (+) a b

I find this code quite inelegant, as it doesn't capture the essence of the vector algebra that's happening (and is also more inefficient since I'm adding and subtracting S[iz]).

One solution would be to implement some kind of vector space typeclass or use a more specific linear algebra library, but since those are such common operations I was wondering if there's a more idiomatic 'straight' solution.

Answer 1

I'd start with a +1 for dfeuer; a more specific library will almost certainly be both cleaner and more efficient.

However, if what you are looking for is a more idiomatic implementation of your centroid function, I like this one:

centroid' :: Simplex -> Int -> Node
centroid' s iz = let t = foldl1 (V.zipWith (+)) (V.drop iz s)
                     n = fromIntegral (V.length t - 1)
                 in V.map (/ n) t

One general comment on your version: it's very easy to create "write-only" Haskell code. There is so much going on in your first line that it is difficult to parse. Your where block is step in the right direction, but I'd go even further to break out the conceptual components.

Also, Hoogle. I didn't know there existed a function drop, but I knew that if it existed, it took an Int and a Vector onto a new Vector. Hoogle doesn't index Vector, but the API for vector is very similar to the API for lists. I searched for "[a] -> Int -> [a]" and "Int -> [a] -> [a]" and found drop.

(Stackage does index Vector, so search for "Int -> Vector a -> Vector a" works there)