1 | function vgg_matrix_test
|
---|
2 | % Test the matrix "flattening" functions from Magnus and Neudecker
|
---|
3 | % vgg_vec
|
---|
4 | % vgg_vech
|
---|
5 | % vgg_commut_matrix
|
---|
6 | % vgg_duplic_matrix
|
---|
7 | % vgg_lmultiply_matrix
|
---|
8 |
|
---|
9 | % Author: awf@robots.ox.ac.uk
|
---|
10 |
|
---|
11 | A = randn(4,5);
|
---|
12 | B = randn(5,2);
|
---|
13 | C = randn(2,4);
|
---|
14 |
|
---|
15 | fprintf('vgg_matrix_test: BEGIN\n');
|
---|
16 | assert('vgg_vec(A*B*C)', 'kron(C'',A)*vgg_vec(B)');
|
---|
17 | assert('vgg_vec(A*B)', 'vgg_lmultiply_matrix(A,size(B))*vgg_vec(B)');
|
---|
18 | assert('vgg_vec(B*C)', 'vgg_rmultiply_matrix(C,size(B))*vgg_vec(B)');
|
---|
19 |
|
---|
20 |
|
---|
21 | % Now try to solve A X + X B = C:
|
---|
22 | X_true = randn(4,4);
|
---|
23 | A = randn(4,4);
|
---|
24 | B = randn(4,4);
|
---|
25 | C = A*X_true + X_true*B;
|
---|
26 |
|
---|
27 | % Re-express in terms of x = vec(X):
|
---|
28 | % fA x + gB x = vec(C)
|
---|
29 | % so
|
---|
30 | % x = (fA + gB) \ vec(C)
|
---|
31 | size_X = [4 4]; % need to know size(X) to convert the matrix multiplys
|
---|
32 | fA = vgg_lmultiply_matrix(A, size_X);
|
---|
33 | gB = vgg_rmultiply_matrix(B, size_X);
|
---|
34 | x = (fA + gB)\vgg_vec(C);
|
---|
35 |
|
---|
36 | % Now check we got it right...
|
---|
37 | assert('x', 'vgg_vec(X_true)')
|
---|
38 |
|
---|
39 |
|
---|
40 |
|
---|
41 | fprintf('vgg_matrix_test: END\n');
|
---|
42 |
|
---|
43 | function assert(sx,sy)
|
---|
44 | x = evalin('caller', sx);
|
---|
45 | y = evalin('caller', sy);
|
---|
46 | err = norm(x-y);
|
---|
47 | if err > 1e-12
|
---|
48 | dbstack
|
---|
49 | fprintf('vgg_matrix_test: FAILED %s = %s\n', sx, sy);
|
---|
50 | keyboard
|
---|
51 | else
|
---|
52 | fprintf('vgg_matrix_test: PASSED %s = %s\n', sx, sy);
|
---|
53 | end
|
---|