[37] | 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
|
---|