Blackbox Examples with Singular Values

This script demonstrates how BBSVDS can be used to compute singular values. To start it directly, use:

 playshow bbdemo_svd1

Construct test-operator
Compute singular values
Reducing memory usage
Resource comparison
Product counts

Construct test-operator

First, we generate a small matrix with singular values that eventually decreases exponentially. To make it a bit challenging, we permute the matrix and premultiply by a dense orthogonal matrix (a scaled Hadamard matrix).

To illustrate the computational effort, we wrap the matrix inside a counting operator.

m=1600; n=500;
sig0=1.02./(1+.02*exp(.065*(0:499)'));
A=spdiags(sig0,0,m,n);
A=A(randperm(m),randperm(n));
A=bbhadamard(m)*blackbox(A)/sqrt(m);
Ac=bbcount(A);

title('Evaluation-tree for test operator');
bbshowtree(Ac,'parent',gca);
axis([-2.5 2.5 -4.5 .5],'off');

Compute singular values

Now compute the singular values with a call to BBSVDS. For reference we time this and save some count statistics.

t=cputime;
sig1=bbsvds(Ac,500);
T1=cputime-t;

cstat1=bbcountstat(Ac,'bbsvds_count');
bbcountreset(Ac);

l1=length(sig1);
semilogy(sig0,'.b'); hold on; grid on;
semilogy(abs(sig0(1:l1)-sig1),'.g');
title('Accuracy of singular values');
xlabel('index');
ylabel('Singular value and absolute error');

500 triplets converged.

Reducing memory usage

Although this is a toy example, we may imagine that we are limited by memory constraints. We can reduce the memory used by modifying the default options.

Normally, BBSVDS will explicitly deflate all triplets that converged previously. This may require a lot of memory and computation time.

We now limit the amount of deflated triplets to 100+20. When more triplets converge, some of them will be "inflated", i.e. they will be overwritten (see also BBDEMO_SVD2).

opt1=bbsvdopt(A,500,false);
opt2=opt1; opt2.dlead=100; opt2.dtail=20;
t=cputime;
sig2=bbsvds(Ac,500,opt2);
T2=cputime-t;

cstat2=bbcountstat(Ac,'bbsvds_count');
bbcountreset(Ac);

l2=length(sig2);
semilogy(abs(sig0(1:l2)-sig2), '.r');

  Total number of converged triplets: 179
Halting: triplets stopped converging.
Increase dlead or rmaxi if you need more triplets.
185 triplets converged.

Resource comparison

The accuracy of the two methods should be comparable. To check how much memory we needed for the computations, we query BBSVDMEM with the options used furing the computations:

mem1=bbsvdmem(Ac,opt1);
mem2=bbsvdmem(Ac,opt2);
fprintf('Converged:         %5i, %5i\n', l1, l2);
fprintf('Computation time:  %5.2f, %5.2f  (s)\n', T1, T2);
fprintf('Memory usage:      %5.2f, %5.2f  (Mb)\n', mem1, mem2);

Converged:           500,   185
Computation time:   6.73, 14.36  (s)
Memory usage:       8.15,  5.25  (Mb)

Product counts

In the first computation, much computational effort went into keeping vectors mutually orthogonal. On the average, this allows the triplets to be resolved with roughly two matrix-vector multiplications, which is the theoretical lower limit (one with the operator and one with the adjoint).

In the second computation the amount of reorthogonalization was reduced. This allows us to use less memory, but needs more multiplications when the singular triplets are inflated after the first 120 (i.e. dlead+dtail) triplets converged. Ultimately the rounding errors limits the number of triplets that will converge.

hold off
plot([cstat1.mm],'.g'); hold on
plot([cstat2.mm],'.r'); grid on
title('Effort spent in computing triplets')
xlabel('Converged triplets')
ylabel('Matrix-vector products')
bbclear Ac;