Owl_optimise.numerical

let numerical_gradient l p x y y' =
  let open MX in
  let h = 0.00001 in
  let f = l y in
  let g = mapi_by_row ~d:(col_num p)
  (fun i v ->
    let fa = f (x $@ (replace_row (v -$ h) p i)) in
    let fb = f (x $@ (replace_row (v +$ h) p i)) in
    (fb -@ fa) /$ (2. *. h)
  ) p in g