1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
use Test::More tests => 13;
use File::Spec;
use lib File::Spec->catfile("..","lib");
use Math::MatrixReal;
do 'funcs.pl';
my $DEBUG2 = 0;
my $string = "[ 1 2 3 ]\n[ 2 2 -1 ]\n[ 1 1 1 ]\n";
my $matrix33 = Math::MatrixReal->new_from_string($string);
print "$matrix33" if $DEBUG;
#
# Tests eigenvalues & eigenvectors computation
#
#
# First the tridiagonal reduction (Householder) on the 3x3
#
my $symm = $matrix33 + ~$matrix33;
ok( not($matrix33->is_symmetric()));
ok( $symm->is_symmetric(), 'A + ~A is symmetric');
print "Matrix 3x3 for eigenvalues & eigenvectors computation:\n$symm"
if $DEBUG;
print "Householder reduction...\n" if $DEBUG;
my ($T, $Q) = $symm->householder();
print "T=\n$T Q=\n$Q" if $DEBUG2;
print "Is Q orthogonal?\n" if $DEBUG;
print ($Q * ~$Q) if $DEBUG2;
ok_matrix_orthogonal($Q);
ok_matrix( $symm, $Q * $T * ~$Q, 'symmetric householder reduction works');
print "Diagonalization of tridiagonal...\n" if $DEBUG;
my ($L1, $V1) = $T->tri_diagonalize($Q);
print "eigenvalues L:\n$L1 eigenvectors:\n$V1" if $DEBUG2;
ok_eigenvectors($symm, $L1, $V1);
# Get first eigenvector
my $aev1 = $V1->column(1);
my $al1 = $L1->element(1,1);
my $ap1_1 = $symm * $aev1; # A * x
my $ap1_2 = $al1 * $aev1; # lambda *x
print "Original computation of A*ev1:\n$ap1_1 Scaled eigenvector:\n$ap1_2"
if $DEBUG2;
ok_matrix( $ap1_1, $ap1_2, 'eigenvectors match');
print "Direct diagonalization...\n" if $DEBUG;
my ($L12, $V12) = $symm->sym_diagonalize();
print "eigenvalues L:\n$L12 eigenvectors:\n$V12" if $DEBUG2;
ok_eigenvectors($symm, $L12, $V12);
ok_matrix_orthogonal($V12);
# Double check the equality
ok_matrix( $L12, $L1);
ok_matrix( $V12, $V1);
#
# Now test the eigenvalues only computations...
#
print "Recomputing: Eigenvalues only.\n 3x3\n" if $DEBUG;
my $altT = $symm->householder_tridiagonal();
ok_matrix( $altT, $T,'householder_tridiagonal works');
my $altL1 = $altT->tri_eigenvalues();
ok_matrix( $altL1, $L1,'tri_eigenvalues works');
my $altL12 = $symm->sym_eigenvalues();
ok_matrix( $altL12, $L12, 'sym_eigenvalues works');
__END__
# Attempt:
# Obtain the eigenvectors when eigenvalues are known
# using inverse iteration.
# We solve (M - lambda * I) * b(k+1) = b(k)
# with b(0) a random unit vector and b(k) is
# normalized at each step.
# This should converge towards the eigenvector,
# but there are problems:
# - for some value, there is not convergence
# (the above system is rather singular, so...)
# - the solution can oscillate between v and -v (?)
# Rodolphe Ortalo, 99/06/14
sub obtain_eigenvector ($$)
{
my ($M, $eigenvalue) = @_;
# Form the linear system A - lamda1 * I
my $inv_it = $M->shadow();
$inv_it->one();
$inv_it->multiply_scalar($inv_it, (-1.0 * $eigenvalue));
$inv_it->add($M, $inv_it);
print "Linear system matrix:\n $inv_it" if $DEBUG2;
# Creates a random vector
my ($rows, $cols) = $inv_it->dim();
my $b = Math::MatrixReal->new($rows, 1);
for (my $i = 1; $i <= $rows; $i++)
{
$b->assign($i, 1, rand());
}
# Normalize it
my $l = $b->length();
$b->multiply_scalar($b, (1.0 / $l));
# Now do LR decomposition for linear system
my $inv_it_LR = $inv_it->decompose_LR();
# Check iterations
my $iter = 0;
my $delta;
do {
my ($dim, $b_base, $base) = $inv_it_LR->solve_LR($b);
# Normalize
my $l = $b_base->length();
$b_base->multiply_scalar($b_base, (-1.0 / $l));
# print "b_base=\n$b_base";
$b->subtract($b_base,$b);
$delta = $b->norm_one();
print "delta=$delta\n";
$b = $b_base;
} while (($delta >= 1e-10) && ($iter++ <= 10));
return $b;
}
#
# Now, try to find one eigenvector again...
# (Using Steffen's functions...:-)
#
my $ev = obtain_eigenvector($symm, $al1);
print "Real ev:\n $aev1 Found ev:\n $ev" if $DEBUG;
ok_matrix(15, $ev, $aev1);
ok_matrix(16, obtain_eigenvector($symm, $L1->element(2,1)),
$V1->column(2));
ok_matrix(17, obtain_eigenvector($symm, $L1->element(3,1)),
$V1->column(3));
|
