Strange incorrect results & different results for identical calls

I was implementing something more complicated with hmatrix & hmatrix-glpk and I ran into issues with the GLPK calls labeling something unexpectedly as Infeasible or Unbounded. I added trace statements to dump the input data (which was just two lists of Doubles), and if I ran this same data on its own in another context (like ghci), it would invariably return Optimal. This continued even to the point where I encoded the input via Data.Binary and then decoded it elsewhere and evaluated it manually there with exactly the same function.

Here is the basic code that causes the issue. I'm sorry it's convoluted - I can't seem to reduce it further without failing to demonstrate the bug, and it's a heavily-pruned version of the algorithm that produced the problem:

module Main where

import Debug.Trace (trace)
import Numeric.LinearAlgebra
import Numeric.LinearProgramming

count = 25
-- Bug is present only with count >= 25
p0 = (2><count) $ repeat 1

-- Bug is not present if I change [Matrix R] and [Double R], which I
-- use only as lists with two elements, to tuples
recurse :: Int -> Int -> [Matrix R] -> [Matrix R]
recurse _ 0 m = m
recurse x iters m = recurse x (iters - 1) d
  where a = f $ map (\q -> flatten $ q ? [x]) m
        d = zipWith (\v q -> accum q (\_ _ -> v) [((x, 1), v)]) a m

f :: [Vector R] -> [Double]
f [p, q] = let obj = zipWith (+) (toList p) (toList q)
               -- Addition above seems necessary for the bug
               constr = Dense $ [ replicate count 1 :==: 1 ]
               bounds = map (:&: (0,1)) [1..count]
               tryA = case simplex (Maximize obj) constr bounds of
                 Optimal  (s, l) -> let r = vector l in [p <.> r, q <.> r]
                 Feasible (s, l) -> let r = vector l in [p <.> r, q <.> r]
                 -- 'Feasible' case must be present above for bug
                 t@_ -> error ("Failed at try A: " ++ show t ++ ", obj=" ++ show obj)
               tryB = case simplex (Maximize obj) constr bounds of
                 Optimal  (s, l) -> let r = vector l in [p <.> r, q <.> r]
                 t@_ -> error ("Failed at try B: " ++ show t ++ ", obj=" ++ show obj)
  in seq tryA tryB

main :: IO ()
main = do
  putStrLn $ show $ recurse 1 100 [p0, p0]

This is also available with a cabal & stack file at https://github.com/Hodapp87/hmatrix-test. I observe the when running stack build && stack exec test:

test: Failed at try B: Unbounded, obj=[2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0]

I do not observe this same behavior with ghci.

Take a look at a few things around line 19: It's solving over 25 variables, it has a constraint that they all sum to 1, they're all bounded between 0 and 1. In the output above, obj is coefficients for the objective function to maximize. It can't be unbounded - it's just a convex sum over a bunch of 2s, its maximum is 2. Evaluating in ghci with the same code and obj value readily shows this:

*Main> let constr = Dense $ [ replicate count 1 :==: 1 ]
*Main> let bounds = map (:&: (0,1)) [1..count]
*Main> let obj=[2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0]
*Main> simplex (Maximize obj) constr bounds
Optimal (2.0,[1.0,-0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0])

Furthermore: in lines 24 and 28, tryA and tryB have identical inputs, yet they produce different results if I force evaluation of both (the error is never produced in evaluating tryA). If I replace seq tryA tryB with just tryA or tryB, the bug does not appear.

So, why are two identical calls producing different results, one of which is incorrect?

haskell-numerics / hmatrix

Strange incorrect results & different results for identical calls #195