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

Revmove reference to gpl2.txt and further shorten copyright and license statement

File size: 2.5 KB
Line 
1% -------------------------------------------------------------------
2% ckbs: Constrained Kalman-Bucy Smoother Program: Copyright (C) 2006
3% Authors: Bradlely Bell:        bradbell at washington dot edu
4%          Gianluigi Pillonetto: giapi at dei dot unipd dot it
5% License: GNU General Public License Version 2
6% -------------------------------------------------------------------
7% $begin sumsq_hes_ok.m$$ $newlinech %$$
8% $spell
9%       ckbs
10%       dg
11%       dh
12%       diff
13%       end end
14%       gradm
15%       gradp
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%
32% $head Source Code$$
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;
82                gradm = ckbs_sumsq_grad(xm, z, g, h, dg, dh, qinv, rinv);
83                %
84                xp       = x;
85                xp(i, k) = xp(i, k) + step;
86                gradp = ckbs_sumsq_grad(xp, z, g, h, dg, dh, qinv, rinv);
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.