Apr
10

Jacobi法求特征值和特征向量[perl]

Categories: Perl | Tags: , | Views: 1,319

 

这个矩阵要求是实对称矩阵，jacobi法的实质就是坐标旋转。对称矩阵和二次型是对应的，通过坐标旋转可以消去交叉项，将原矩阵化成只剩对角元素的三角阵（其他元素为0），这些对角元素就是矩阵的特征值。

可以证明，jacobi方法是收敛的。其缺点是对于稀疏矩阵旋转后难保持其稀疏性。

具体的理论算法请看这里

 

 

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
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229

#
# Panwen Wang, April 8th, 2009
# pwwang AT pwwang.com
# get eigenvalues and eigenvectors of
# a real symmetric matrix like
#  1  -1  0 
# -1   1  2
#  0   2  1
#
# Algorithm:
# 1. get a nonzero and non-diagonal element a[i,j] of matrix A, usually
#    the maximun absolute value of non-diagonal elements of A (max(A));  
# 2. give sin(phi) and cos(phi) by 
#    ( a[j,j]-a[i,i] )*sin(2*phi) + 2*a[i,j]*cos(2*phi) = 0 ;
#    => tan(2*phi) = 2*a[i,j]/(a[i,i]-a[j,j]) 
#    => phi = arctan(2*a[i,j]/(a[i,i]-a[j,j]))/2 
# 3. get elements of a new matrix A1(a1[i,j]) by
#    a1[i,i] = a[i,i]*cos^2(phi) + a[j,j]*sin^2(phi) + 2*a[i,j]*cos(phi)*sin(phi)
#    a1[j,j] = a[i,i]*sin^2(phi) + a[j,j]*cos^2(phi) - 2*a[i,j]*cos(phi)*sin(phi)
#    a1[i,l] = a1[l,i] = a[i,l]*cos(phi) + a[j,l]*sin(phi)       ( l!=i,j )
#    a1[j,l] = a1[l,j] = -a[i,l]*sin(phi) + a[j,l]*cos(phi)      ( l!=i,j )
#    a1[l,m] = a1[m,l] = a[m,l]                                  ( m,l != i,j )
#    a1[i,j] = a1[j,i] = { (a[j,j]-a[i,i])*sin(2*phi) }/2 + a[i,j]*(cos^2(phi) - sin^2(phi))
# 4. let A1 be the substitution of A, repeat step 1,2,3 and get A2, and A3,A4,...,An can be 
#    obtained by the same way. Calculation ceases if max(An) is less than the given threshold.
#
package Jacobi;
use strict;
 
sub new{
    my $class     = shift;
    my $self    = {
        MATRIX        => [],
        EIGENVALUE  => [],
        EIGENVECTOR => [],
        PRECISION       => 1e-5,
        DIMENSION   => undef,
    };
    bless $self,$class;
    return $self;
}
 
# load matrix
sub setMatrix{
    my $self = shift;
    $self->{MATRIX} = [ @_ ];
    $self->{DIMENSION} = scalar(@{$self->{MATRIX}});
    $self->jacobi;
}
 
# load matrix from an array
sub setMatrixFromArray{
    my $self = shift;
    my @matrix = ();
    my @array = @_;
    my $s = sqrt(@array);
    if($s =~ /^[1-9]d*(.0*)?$/){
    # tell if $s is an integer.
    # if it is, the array can be converted to a matrix
    # and $s is the dimension of the matrix
    # if not, exit
        $self->{DIMENSION} = $s;
        for(my $i=0; $i<$s; $i++){
            my @row = ();
            for(my $j=0; $j<$s; $j++){
                push @row, $array[$i*$s+$j];
            }
            push @matrix, [ @row ];
        }
    } else {
        die "setMatrixFromArray: Array cannot be converted to matrix.n";
    }
    $self->{MATRIX} = [ @matrix ];
    $self->jacobi;
}
 
# set precision request
sub setPrecision{
    my $self = shift;
    $self->{PRECISION} = shift;
}
 
# eigenvalues
sub eigenvalues{
    my $self = shift;
    return @{$self->{EIGENVALUE}};
}
 
# eigenvectors
sub eigenvectors{
    my $self = shift;
    return @{$self->{EIGENVECTOR}};
}
 
# get max element fabs and its subscripts
sub max{
    my $self = shift;
    my $matrix = [ @_ ];
    my $max = abs($matrix->[0][1]);
    my ($i, $j);
    my @ijmax = (0,1,$max); 
    for($i=0; $i<$self->{DIMENSION}; $i++){ 
    # ensure $i is less than $j
        for($j=$i; $j<$self->{DIMENSION}; $j++){
            if( $i!=$j ){
                if( abs($matrix->[$i][$j]) > $max){
                    $max = abs($matrix->[$i][$j]);
                    @ijmax = ($i, $j, $max);
                } # if
            } # if
        } # for $j
    } # for $i
    return @ijmax;
}
 
# step 3, get a new matrix by a give phi value
sub newMatrix{
    my $self = shift;
    my ($i, $j, $phi, @mat) = @_;
    my $sp = sin($phi);
    my $cp = cos($phi);
    my @mat1 = ();
 
    for(my $ii=0; $ii<$self->{DIMENSION}; $ii++){
        for(my $jj=$ii; $jj<$self->{DIMENSION}; $jj++){
            if( $ii==$i ){      # row $i
                if( $jj==$i ){      # colomn $i
                    $mat1[$ii][$jj] = $mat[$i][$i]*$cp*$cp + $mat[$j][$j]*$sp*$sp + 2*$mat[$i][$j]*$cp*$sp; 
                } elsif( $jj==$j ){ # colomn $j
                    $mat1[$ii][$jj] = ($mat[$j][$j]-$mat[$i][$i])*$sp*$cp + $mat[$i][$j]*($cp*$cp-$sp*$sp); 
                } else {            # colomn $l ( l!=i,j)
                    $mat1[$ii][$jj] = $mat[$i][$jj]*$cp + $mat[$j][$jj]*$sp;
                }
            } elsif ( $ii==$j ){# row $j
                if( $jj==$i ){      # colomn $i
                    $mat1[$ii][$jj] = ($mat[$j][$j]-$mat[$i][$i])*$sp*$cp + $mat[$i][$j]*($cp*$cp-$sp*$sp);
                } elsif( $jj==$j ){ # colomn $j
                    $mat1[$ii][$jj] = $mat[$i][$i]*$sp*$sp + $mat[$j][$j]*$cp*$cp - 2*$mat[$i][$j]*$cp*$sp;
                } else {            # colomn $l ( l!=i,j )
                    $mat1[$ii][$jj] = $mat[$j][$jj]*$cp - $mat[$i][$jj]*$sp;
                }
            } else {            # row $l ( l!=i,j )
                if( $jj==$i ){      # colomn $i
                    $mat1[$ii][$jj] = $mat[$i][$ii]*$cp + $mat[$j][$ii]*$sp;
                } elsif( $jj==$j ){ # colomn $j
                    $mat1[$ii][$jj] = $mat[$j][$ii]*$cp - $mat[$i][$ii]*$sp;
                } else {            # colomn $l ( l!=i,j )
                    $mat1[$ii][$jj] = $mat[$ii][$jj];
                }
            }
            $mat1[$jj][$ii] = $mat1[$ii][$jj];
        }
    }    
    return @mat1;
}
 
# calculating
sub jacobi{
    my $self = shift;
    my @matrix = @{$self->{MATRIX}};  # initial matrix
    my $phi;
    my @tempVectors = ();
    my @vectors = ();
    my @cmat = ();
    for(my $x=0; $x<$self->{DIMENSION}; $x++){
        for(my $y=0; $y<$self->{DIMENSION}; $y++){
            if($x==$y) { $tempVectors[$x][$y] = 1.0; } 
            else { $tempVectors[$x][$y] = 0.0; }
        }
    }
 
    for( ;; ){ # step 4
        my ($i, $j, $max) =  $self->max(@matrix);  # step 1
        last if( $max < $self->{PRECISION} ); 
        # if the maximum value of the matrix LE(Less than or Equal to) the limit, break
        $phi = atan2(2*$matrix[$i][$j], $matrix[$i][$i] - $matrix[$j][$j]) / 2;  # step 2
        @matrix = $self->newMatrix($i,$j,$phi,@matrix); # step 3
        for(my $x=0; $x<$self->{DIMENSION}; $x++){
            for(my $y=0; $y<$self->{DIMENSION}; $y++){
                if($x==$y) { $cmat[$x][$y] = 1.0; }
                else { $cmat[$x][$y] = 0.0; }
                $vectors[$x][$y] = 0.0;
            }
        }    
        $cmat[$i][$i] = cos($phi);
        $cmat[$j][$j] = cos($phi);
        $cmat[$i][$j] = -sin($phi);
        $cmat[$j][$i] = sin($phi);    
        for(my $x=0; $x<$self->{DIMENSION}; $x++){ # for eigenvectors
            for(my $y=0; $y<$self->{DIMENSION}; $y++){
                for(my $z=0; $z<$self->{DIMENSION}; $z++){
                    $vectors[$x][$y] = $vectors[$x][$y] + $tempVectors[$x][$z] * $cmat[$z][$y];
                }
                #print $vectors[$x][$y],"-",$tempVectors[$x][$y]," ";
            }
        }
        for(my $x=0; $x<$self->{DIMENSION}; $x++){
            @{$tempVectors[$x]} = @{$vectors[$x]};
        }
        #@tempVectors = @vectors;
    }
    # the digonal elements are the eigenvalues
    for(my $x=0; $x<$self->{DIMENSION}; $x++){
        push @{$self->{EIGENVALUE}}, $matrix[$x][$x];
    }
    @{$self->{EIGENVECTOR}} = ();
    # perl does not have a conception of multiple array, thus it can not give a colomn
    # by a simple variable, so we transpose @vectors, let row be the eigenvector of the
    # corresponding eigenvalue.
    for(my $x=0; $x<$self->{DIMENSION}; $x++){
        for(my $y=0; $y<$self->{DIMENSION}; $y++){
            $self->{EIGENVECTOR}[$x][$y] = $vectors[$y][$x];
        }
    }
}
 
sub printEig{
    my $self = shift;
    for(my $i=0; $i<$self->{DIMENSION}; $i++){
        print "Eigenvalue [$i]: ", $self->{EIGENVALUE}[$i], "n";
        print "Eigenvector[$i]: [ ";
        for(my $j=0; $j<$self->{DIMENSION}; $j++){
            printf "%9.6f, ",$self->{EIGENVECTOR}[$i][$j];
        }
        print " ] n";
    }
}
 
1;

 

 

调用方法：

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24

#!/usr/bin/perl
 
use strict;
use Jacobi;  # package name is Jacobi.pm
 
my @a = (
    [2,-1,0],
    [-1,2,-1],
    [0,-1,2]
);
 
#my @a = (
#    [10,     7,     8,     7],
#    [ 7,     5,     6,     5],
#    [ 8,     6,    10,     9],
#    [ 7,     5,     9,    10]
#);
 
my $jac = Jacobi->new; 
$jac->setMatrix(@a);
my @eigenvalues = $jac->eigenvales;
my @eigenvectors = $jac->eigenvectors; # 2-D array
# each ROW corresponds to the corresponding eigenvalue.
$jac->printEig;

 

这篇文章来自 迷途知返(PWWANG.COM), 转载请注明出处。 版权说明

Leave a comment | Trackback
  1. Rock
    April 11th, 2009 at 08:40
    Reply | Quote | #1

    呵呵,这几天我也在研习3D数学了,公司有个引擎似乎要完善了…

;) :| :x :twisted: :roll: :oops: :o :mrgreen: :lol: :idea: :evil: :cry: :arrow: :P :D :?: :? :) :( :!: 8O 8)

你可以使用@somebody:开头, 来邮件通知somebody你回复了他的留言(用户名区分大小写).