# source:projects/ckbs/trunk/test/sumsq_hes_ok.m@95

Last change on this file since 95 was 35, checked in by bradbell, 11 years ago

File size: 2.5 KB
Line
1% -------------------------------------------------------------------
2% ckbs: Constrained Kalman-Bucy Smoother Program: Copyright (C) 2006
4%          Gianluigi Pillonetto: giapi at dei dot unipd dot it
6% -------------------------------------------------------------------
7% \$begin sumsq_hes_ok.m\$\$ \$newlinech %\$\$
8% \$spell
9%       ckbs
10%       dg
11%       dh
12%       diff
13%       end end
16%       hes
17%       qinv
18%       rinv
19%       sumsq
20%       tmp
21%       xm
22%       xp
23% \$\$
24%
25% \$section ckbs_sumsq_hes Example and Test\$\$
26%
27% \$index ckbs_sumsq_hes, example and test\$\$
28% \$index sumsq_hes, example and test\$\$
29% \$index example, sumsq_hes\$\$
30% \$index test, sumsq_hes\$\$
31%
33% \$newlinech \$\$ \$codep
34function [ok] = sumsq_hes_ok()
35ok = true;
36% --------------------------------------------------------
37% You can change these parameters
38m    = 1;   % number of measurements per time point
39n    = 2;   % number of state vector components per time point
40N    = 3;   % number of time points
41% ---------------------------------------------------------
42%  Define the problem
43rand('seed', 123);
44dg   = zeros(n, n, N);
45dh   = zeros(m, n, N);
46qinv = zeros(n, n, N);
47rinv = zeros(m, m, N);
48for k = 1 : N
49        dh(:, :, k)   = rand(m, n);
50        dg(:, :, k)   = rand(n, n);
51        tmp           = rand(m, m);
52        rinv(:, :, k) = (tmp + tmp') / 2 + 2 * eye(m);
53        tmp           = rand(n, n);
54        qinv(:, :, k) = (tmp + tmp') / 2 + 2 * eye(n);
55end
56% ---------------------------------------------------------
57% Compute the Hessian using ckbs_sumsq_hes
58[D, A] = ckbs_sumsq_hes(dg, dh, qinv, rinv);
59% ---------------------------------------------------------
60H    = zeros(n * N , n * N );
61for k = 1 : N
62        index           = (k - 1) * n + (1 : n);
63        H(index, index) = D(:, :, k);
64        if k > 1
65                H(index - n, index) = A(:, :, k)';
66                H(index, index - n) = A(:, :, k);
67        end
68end
69%
70% Use finite differences to check Hessian
71x    = rand(n, N);
72z    = rand(m, N);
73h    = rand(m, N);
74g    = rand(n, N);
75%
76step   = 1;
77for k = 1 : N
78        for i = 1 : n
79                % Check second partial w.r.t x(i, k)
80                xm       = x;
81                xm(i, k) = xm(i, k) - step;
83                %
84                xp       = x;
85                xp(i, k) = xp(i, k) + step;
87                %
88                check  = (gradp - gradm) / ( 2 * step );
89                for k1 = 1 : N
90                        for i1 = 1 : n
91                                value = H(i + (k-1)*n, i1 + (k1-1)*n);
92                                diff  = check(i1, k1) - value;
93                                ok     = ok & ( abs(diff) < 1e-10 );
94                        end
95                end
96        end
97end
98return
99end
100% \$\$ \$newlinech %\$\$
101% \$end
Note: See TracBrowser for help on using the repository browser.