[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
normest for octave
From: |
David Bateman |
Subject: |
normest for octave |
Date: |
Mon, 04 Dec 2006 16:58:15 +0100 |
User-agent: |
Thunderbird 1.5.0.7 (X11/20060921) |
I've been trying to get the last couple of missing functions for sparse
functionality written, ignoring the functions that are obsolete such as
symmmd and apart from the iterative solvers. These missing functions
are normest, condest and svds.. I'll wait for svds till eigs is in
octave, as svds basically uses eigs on [sparse(size(A,1),size(A,1)), A;
A', sparse(size(A,2),size(A,2))] together with an estimate of the 2-norm
of the matrix.
So the first of the three missing functions I think needs to be
addressed in normest. The power series for finding the largest
eigenvector and eigenvalue can be written as V1 = Lambda * A * V0; and
for N iterations might be written as
V = randn(size(A,2), 1); V = V / norm(V); for i=1:N, t = A * V; Lambda =
norm(t); V = t / Lambda; end
So to get the largest singular value I can use the same trick as for
svds and use
B = [sparse(size(A,1),size(A,1)), A; A', sparse(size(A,2),size(A,2))] ;
V = randn(size(B,2), 1); V = V / norm(V); for i=1:N, t = B * V; Lambda =
norm(t); V = t / Lambda; end
where Lambda then returns an estimate of the 2-norm of the matrix A.
The attached function, with a few extras, does exactly the above
operation. However, I've found that it takes more iterations than the
equivalent matlab function to converge (about twice as many). Does
anyone know how I can accelerate the attached functions convergence? Can
I do it without explicitly forming the matrix B?
Regards
David
--
David Bateman address@hidden
Motorola Labs - Paris +33 1 69 35 48 04 (Ph)
Parc Les Algorithmes, Commune de St Aubin +33 6 72 01 06 33 (Mob)
91193 Gif-Sur-Yvette FRANCE +33 1 69 35 77 01 (Fax)
The information contained in this communication has been classified as:
[x] General Business Information
[ ] Motorola Internal Use Only
[ ] Motorola Confidential Proprietary
function [e1, c] = mynormest(A, tol)
if (nargin < 2)
tol = 1e-6;
end
if (tol < eps)
tol = eps
endif
if (ndims(A) != 2)
error("A must be a matrix");
endif
[m, n] = size(A);
B = [sparse(m,m), A; A', sparse(n,n)];
c = 0;
x = sum(abs(B))';
e0 = 0;
e1 = norm(x);
x = x / e1;
while (abs(e1 - e0) > tol)
e0 = e1;
x = B * x;
e1 = norm(x);
x = x / e1;
c++;
endwhile
endfunction
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- normest for octave,
David Bateman <=