On Thu, 2 Mar 2006, John Randall wrote:
> Here is a slightly more comprehensible form of the power method code I
> posted yesterday. It will correctly calculate the principal eigenvectors
> in Coyle's paper.
>
> mp=:+/ . *
> normalize=:%(>./@:|)
> iterate=:normalize@:mp
> init=:[: ? # # 0:
> power=:13 : 'y.&iterate^:_ init y.'
>
> rod=. >(1,(%3),%9);(3,1,%3);9 3 1 NB. "Replication of detail", p. 7
>
> (%+/) power rod
> 0.0769231 0.230769 0.692308
>
> 0.077 0.231 0.692 NB. Coyle
The iteration might in general be accelerated if you square the matrix
at each iteration (and adjust to avoid possible eventual numeric
overflow), since M, M2=. M mp M, M4=. M2 mp M2... all have the
same eigenvectors. This is useful if the 2nd largest eigenvalue
is close to the principal eigenvalue.
Quick hack...
%----------------------------- J script -----------------------------
mp=: +/ . *
init=: [: ? # # 0:
power2=: 3 : 0
v0=. init y.
while. v0 ~: v1=. y. mp v0 do.
y.=. (mp~ y.) % l2=. +/@:(* +) v0=. v1
v0=. v0 % %: l2
end.
)
rod=: >(1,(%3),%9);(3,1,%3);9 3 1 NB. "Replication of detail", p. 7
%----------------------------- J example -----------------------------
power2 rod
0.104828 0.314485 0.943456
%----------------------------- end J -----------------------------
Note the normalisation is different (eigenvector has length 1)
since this is cannibalised from one of my antediluvian scripts
[treating the availability of LAPACK in J as the deluge :-)]
-- Ewart Shaw
J.E.H.Shaw [Ewart Shaw] [EMAIL PROTECTED] TEL: +44 2476 523069
Department of Statistics, University of Warwick, Coventry CV4 7AL, UK
http://www.warwick.ac.uk/statsdept http://www.ewartshaw.co.uk
3 ((4&({*.(=+/))++/=3:)@([:,/0&,^:(i.3)@|:"2^:2))&.>@]^:([EMAIL PROTECTED])
<#:3 6 2
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
