Source file `owl_numdiff_generic.ml`

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
# 1 "src/base/optimise/owl_numdiff_generic.ml"
(*
 * OWL - an OCaml numerical library for scientific computing
 * Copyright (c) 2016-2018 Liang Wang <liang.wang@cl.cam.ac.uk>
 *)

open Owl_types


(* Functor of making numerical differentiation module of different precisions *)

module Make
  (A : Ndarray_Numdiff)
  = struct

  type arr = A.arr
  type elt = A.elt


  (* global epsilon value used in numerical differentiation *)
  let _eps = 0.00001

  let _ep1 = 1.0 /. _eps

  let _ep2 = 0.5 /. _eps


  (* derivative of f : scalar -> scalar *)
  let diff f x = (f (x +. _eps) -. f (x -. _eps)) *. _ep2


  (* derivative of f : scalar -> scalar, return both function value and derivative *)
  let diff' f x = f x, diff f x


  (* second order derivative of f : scalar -> scalar *)
  let diff2 f x = (f (x +. _eps) +. f (x -. _eps) -. (2. *. f x)) /. (_eps *. _eps)


  (* second order derivative of f : float -> float, return both function value and derivative *)
  let diff2' f x = f x, diff2 f x


  (* gradient of f : vector -> scalar, return both function value and gradient *)
  let grad' f x =
    let n = A.numel x in
    let g = A.create [|n|] (f x) in
    let gg = A.mapi (fun i xi ->
      let x' = A.copy x in
      A.set x' [|i|] (xi +. _eps);
      f x'
    ) x
    in
    g, A.((gg - g) *$ _ep1)


  (* gradient of f : vector -> scalar *)
  let grad f x = grad' f x |> snd


  (* transposed jacobian of f : vector -> vector, return both function value and jacobian *)
  let jacobianT' f x =
    let y = f x in
    let m, n = A.numel x, A.numel y in
    let j = A.tile y [|m; 1|] in
    let jj = A.copy j in

    for i = 0 to m - 1 do
      let x' = A.copy x in
      A.set x' [|i|] ((A.get x [|i|]) +. _eps);
      let y' = A.reshape (f x') [|1; n|] in
      A.set_slice [[i];[]] jj y'
    done;

    y, A.((jj - j) *$ _ep1)


  (* transposed jacobian of f : vector -> vector *)
  let jacobianT f x = jacobianT' f x |> snd


  (* jacobian of f : vector -> vector, return both function value and jacobian *)
  let jacobian' f x =
    let y, j = jacobianT' f x in
    y, A.transpose ~axis:[|1;0|] j


  (* jacobian of f : vector -> vector *)
  let jacobian f x = jacobian' f x |> snd


end


(* ends here *)