octave-maintainers
[Top][All Lists]
Advanced

[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

reply via email to

[Prev in Thread] Current Thread [Next in Thread]