hessian/
hessian.rs

1// SPDX-License-Identifier: EPL-2.0 OR GPL-2.0-or-later
2// SPDX-FileCopyrightText: Bradley M. Bell <bradbell@seanet.com>
3// SPDX-FileContributor: 2025 Bradley M. Bell
4//
5// Example Hessian
6//
7use rustad::{
8    AD,
9    ad_from_value,
10    ad_from_vector,
11    NumVec,
12    start_recording,
13    stop_recording,
14};
15//
16// example_hessian
17// Simple case where V = f32
18fn example_hessian () {
19    //
20    type V     = f32;
21    let nx     = 3;
22    let trace  = false;
23    //
24    // x
25    let x  : Vec<V> = vec![ 2.0 as V; nx ];
26    //
27    // ax
28    let ax       = start_recording(x.clone());
29    //
30    // asum
31    let mut asum : AD<V>  = ad_from_value(  0.0 as V );
32    for j in 0 .. nx {
33        let cubed  = &( &ax[j] * &ax[j] ) * &ax[j];
34        asum      += &cubed;
35    }
36    //
37    // f
38    // f(x) = x[0] * x[0] * x[0] + ... + x[nx-1] * x[nx-1] * x[nx-1]
39    let ay = vec![ asum ];
40    let f  = stop_recording(ay);
41    //
42    // av
43    let ax                     = start_recording(x);
44    let mut av  : Vec< AD<V> > = Vec::new();
45    f.forward_zero_ad(&mut av, ax, trace);
46    //
47    // g
48    // g(x) = df/dx = [ 3 * x[0] * x[0], ..., 3 * x[nx-1] * x[nx-1] ]
49    let dy  : Vec<V>  = vec![ 1.0 as V ];
50    let ady           = ad_from_vector(dy);
51    let adx           = f.reverse_one_ad(&av, ady, trace);
52    let g             = stop_recording(adx);
53    //
54    // x
55    // x[j] = j+2
56    let mut x  : Vec<V> = Vec::new();
57    for j in 0 .. nx {
58        x.push( (j+2) as V );
59    }
60    //
61    // v, y
62    let mut v  : Vec<V> = Vec::new();
63    let y               = g.forward_zero_value(&mut v, x, trace);
64    for j in 0 .. nx {
65        let check  = 3 * (j+2) * (j+2);
66        assert_eq!( y[j], check as V );
67    }
68    //
69    // dy
70    // dy[i] = partial g[i] w.r.t x[j] = 6 * x[j]
71    for j in 0 .. nx {
72        let mut dx : Vec<V> = vec![ 0.0 as V; nx ];
73        dx[j]               = 1.0 as V;
74        let dy              = g.forward_one_value(&v, dx, trace);
75        for i in 0 .. nx {
76            if i == j {
77                let check  = 6 * (j+2);
78                assert_eq!( dy[i], check as V );
79            } else {
80                assert_eq!( dy[i],  0.0 as V );
81            }
82        }
83    }
84}
85//
86// example_numvec_hessian
87// Same function where V = NumVec<f64>
88fn example_numvec_hessian () {
89    //
90    type F     = f64;
91    type V     = NumVec<F>;
92    let nx     = 3;
93    let trace  = false;
94    //
95    // x
96    let mut x  : Vec<V> = Vec::new();
97    for _j in 0 .. nx {
98        x.push( NumVec::new( vec![ 2.0 as F ] ) );
99    }
100    //
101    // ax
102    let ax       = start_recording(x.clone());
103    //
104    // asum
105    let mut asum : AD<V>  = ad_from_value(  NumVec::from( 0.0 as F ) );
106    for j in 0 .. nx {
107        let cubed  = &( &ax[j] * &ax[j] ) * &ax[j];
108        asum      += &cubed;
109    }
110    //
111    // f
112    // f(x) = x[0] * x[0] * x[0] + ... + x[nx-1] * x[nx-1] * x[nx-1]
113    let ay = vec![ asum ];
114    let f  = stop_recording(ay);
115    //
116    // av
117    let ax                     = start_recording(x);
118    let mut av  : Vec< AD<V> > = Vec::new();
119    f.forward_zero_ad(&mut av, ax, trace);
120    //
121    // g
122    // g(x) = df/dx = [ 3 * x[0] * x[0], ..., 3 * x[nx-1] * x[nx-1] ]
123    let dy  : Vec<V>  = vec![ NumVec::from( 1.0 as F ) ];
124    let ady           = ad_from_vector(dy);
125    let adx           = f.reverse_one_ad(&av, ady, trace);
126    let g             = stop_recording(adx);
127    //
128    // x
129    // x[j] = [ j+1, j+2 ]
130    let mut x  : Vec<V> = Vec::new();
131    for j in 0 .. nx {
132        x.push( NumVec::new( vec![ (j+1) as F, (j+2) as F ] ) );
133    }
134    //
135    // v, y
136    let mut v  : Vec<V> = Vec::new();
137    let y               = g.forward_zero_value(&mut v, x, trace);
138    for j in 0 .. nx {
139        //
140        let check  = 3 * (j+1) * (j+1);
141        assert_eq!( y[j].get(0), check as F );
142        //
143        let check  = 3 * (j+2) * (j+2);
144        assert_eq!( y[j].get(1), check as F );
145    }
146    //
147    // dy
148    // dy[i] = partial g[i] w.r.t x[j] = 6 * x[j]
149    for j in 0 .. nx {
150        let mut dx : Vec<V> = vec![ NumVec::from( 0.0 as F ); nx ];
151        dx[j]               = NumVec::from( 1.0 as F );
152        let dy              = g.forward_one_value(&v, dx, trace);
153        for i in 0 .. nx {
154            if i == j {
155                //
156                let check  = 6 * (j+1);
157                assert_eq!( dy[i].get(0), check as F );
158                //
159                let check  = 6 * (j+2);
160                assert_eq!( dy[i].get(1), check as F );
161            } else {
162                for k in 0 .. dy[i].len() {
163                    assert_eq!( dy[i].get(k) ,  0.0 as F );
164                }
165            }
166        }
167    }
168}
169fn main() {
170    example_hessian();
171    example_numvec_hessian();
172}
hessian/hessian.rs

hessian/
hessian.rs
