function [x y] = snbin(A,nmv,m,n) % Stochastic matrix-free binormalization for nonsymmetric real A. % [x y] = snbin(A,nmv,m,n) % A is a matrix or function handle. If it is a function handle, then % v = A(x) returns A*x and v = A(x,'trans') returns A'*x. % nmv is the number of forward and transpose matrix-vector product pairs to % perform. % m,n is the size of the matrix. It is necessary to specify these only if % A is a function handle. % diag(x) A diag(y) is approximately binormalized. % Jan 2010. Algorithm and code by Andrew M. Bradley (ambrad@stanford.edu). % Aug 2010. New omega schedule. op = isa(A,'function_handle'); if(~op) [m n] = size(A); end r = ones(m,1); c = ones(n,1); for(k = 1:nmv) % omega^k alpha = (k - 1)/nmv; omega = (1 - alpha)*1/2 + alpha*1/nmv; % rows s = randn(n,1)./sqrt(c); if(op) y = A(s); else y = A*s; end r = (1-omega)*r/sum(r) + omega*y.^2/sum(y.^2); % columns s = randn(m,1)./sqrt(r); if(op) y = A(s,'trans'); else y = (s'*A)'; end c = (1-omega)*c/sum(c) + omega*y.^2/sum(y.^2); end x = 1./sqrt(r); y = 1./sqrt(c);