Skip to content

Commit

Permalink
'SVD' variant for spy in @hodlr and @hss.
Browse files Browse the repository at this point in the history
  • Loading branch information
@robol
robol committed Oct 2, 2023
1 parent d331301 commit c0e4dbe
Show file tree
Hide file tree
Showing 2 changed files with 93 additions and 12 deletions.
55 changes: 48 additions & 7 deletions @hodlr/spy.m
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,15 @@
function spy(H)
function spy(H, variant)
%SPY Inspect the rank structure of H.
%
% SPY(H) draw a plot displaying the low-rank structure in the off-diagonal
% blocks of H.
%
% SPY(H, 'svd') also plots bars representing the singular values of the
% low-rank blocks, to indicate the decay present in those blocks.

if ~exist('variant', 'var')
variant = 'standard';
end

m = size(H, 1);
n = size(H, 2);
Expand All @@ -16,15 +23,16 @@ function spy(H)
axis([ 0 n 0 m ]);
hold on;

spy_draw_block(H, 0, 0, hodlrrank(H));
spy_draw_block(H, 0, 0, hodlrrank(H), variant);

hold off;

end

function spy_draw_block(H, xoffset, yoffset, rk)
function spy_draw_block(H, xoffset, yoffset, rk, variant)
diag_color = [ 0.1, 0.3, 0.95 ];
offdiag_color = [ 0.95, 0.95, 1.0 ];
svd_color = [1.0, 0.3, 0.4];

if is_leafnode(H)
mm = size(H, 1);
Expand All @@ -35,10 +43,10 @@ function spy_draw_block(H, xoffset, yoffset, rk)

fill(xb, yb, diag_color);
else
spy_draw_block(H.A11, xoffset, yoffset, rk);
spy_draw_block(H.A11, xoffset, yoffset, rk, variant);
mm = size(H.A11, 1);
nn = size(H.A11, 2);
spy_draw_block(H.A22, xoffset + nn, yoffset + mm, rk);
spy_draw_block(H.A22, xoffset + nn, yoffset + mm, rk, variant);

mm2 = size(H.A22, 1);
nn2 = size(H.A22, 2);
Expand All @@ -48,21 +56,54 @@ function spy_draw_block(H, xoffset, yoffset, rk)
xoffset + nn + nn2, xoffset + nn ];
yb = [yoffset, yoffset, yoffset + mm, yoffset + mm ];
thisrk = size(H.U12, 2);

fill(xb, yb, offdiag_color);

if strcmp(variant, 'svd')
% Draw the SVD plot of the block in the background
[~, R] = qr(H.U12, 0); [~, S] = qr(H.V12, 0);
s = log(svd(R * S'));
nsvd = length(s);
for j = 1 : nsvd
sjm = (s(1) - s(j)) / (s(1) - s(end)) * mm;
xbj = [xoffset + nn + nn2 * (j-1)/nsvd, xoffset + nn + nn2 * j / nsvd, ...
xoffset + nn + nn2 * j / nsvd, xoffset + nn + nn2 * (j-1)/nsvd ];
ybj = [yoffset + sjm, yoffset + sjm, yoffset + mm, yoffset + mm ];
fill(xbj, ybj, svd_color, 'LineStyle','none');
end
end

text(xoffset + nn + nn2 / 2, yoffset + mm / 2, ...
sprintf('%d', thisrk), ...
'HorizontalAlignment', 'center');

% Low rank block (2, 1)
xb = [xoffset, xoffset + nn, xoffset + nn, xoffset ];
yb = [yoffset + mm, yoffset + mm, ...
yoffset + mm + mm2, yoffset + mm + mm2 ];
thisrk = size(H.U21, 2);

fill(xb, yb, offdiag_color);

if strcmp(variant, 'svd')
% Draw the SVD plot of the block in the background
[~, R] = qr(H.U21, 0); [~, S] = qr(H.V21, 0);
s = log(svd(R * S'));
nsvd = length(s);
for j = 1 : nsvd
sjm = (s(1) - s(j)) / (s(1) - s(end)) * mm2;
xbj = [xoffset + nn * (j-1)/nsvd, xoffset + nn * j / nsvd, ...
xoffset + nn * j / nsvd, xoffset + nn * (j-1)/nsvd ];
ybj = [yoffset + mm + sjm, yoffset + mm + sjm, yoffset + mm + mm2, yoffset + mm + mm2 ];
fill(xbj, ybj, svd_color, 'LineStyle','none');
end
end

text(xoffset + nn / 2, yoffset + mm + mm2 / 2, ...
sprintf('%d', thisrk), ...
'HorizontalAlignment', 'center');
end
end

function spy_draw_svd_block()
end
50 changes: 45 additions & 5 deletions @hss/spy.m
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,19 @@
function spy(H)
function spy(H, variant)
%SPY Inspect the rank structure of H.
%
% SPY(H) draw a plot displaying the low-rank structure in the off-diagonal
% blocks of H.
%
% SPY(H, 'svd') also plots bars representing the singular values of the
% low-rank blocks, to indicate the decay present in those blocks.

if ~exist('variant', 'var')
variant = 'standard';
end

if strcmp('variant', 'svd')
H = hss_proper(H);
end

m = size(H, 1);
n = size(H, 2);
Expand All @@ -16,15 +27,16 @@ function spy(H)
axis([ 0 n 0 m ]);
hold on;

spy_draw_block(H, 0, 0, hodlrrank(H));
spy_draw_block(H, 0, 0, hssrank(H), variant);

hold off;

end

function spy_draw_block(H, xoffset, yoffset, rk)
function spy_draw_block(H, xoffset, yoffset, rk, variant)
diag_color = [ 0.1, 0.3, 0.95 ];
offdiag_color = [ 0.95, 0.95, 1.0 ];
svd_color = [1.0, 0.3, 0.4];

if H.leafnode == 1
mm = size(H, 1);
Expand All @@ -35,10 +47,10 @@ function spy_draw_block(H, xoffset, yoffset, rk)

fill(xb, yb, diag_color);
else
spy_draw_block(H.A11, xoffset, yoffset, rk);
spy_draw_block(H.A11, xoffset, yoffset, rk, variant);
mm = size(H.A11, 1);
nn = size(H.A11, 2);
spy_draw_block(H.A22, xoffset + nn, yoffset + mm, rk);
spy_draw_block(H.A22, xoffset + nn, yoffset + mm, rk, variant);

mm2 = size(H.A22, 1);
nn2 = size(H.A22, 2);
Expand All @@ -50,6 +62,20 @@ function spy_draw_block(H, xoffset, yoffset, rk)
thisrk = rank(H.B12);

fill(xb, yb, offdiag_color);

if strcmp(variant, 'svd')
% Draw the SVD plot of the block in the background
s = log(svd(H.B12));
nsvd = length(s);
for j = 1 : nsvd
sjm = (s(1) - s(j)) / (s(1) - s(end)) * mm;
xbj = [xoffset + nn + nn2 * (j-1)/nsvd, xoffset + nn + nn2 * j / nsvd, ...
xoffset + nn + nn2 * j / nsvd, xoffset + nn + nn2 * (j-1)/nsvd ];
ybj = [yoffset + sjm, yoffset + sjm, yoffset + mm, yoffset + mm ];
fill(xbj, ybj, svd_color, 'LineStyle','none');
end
end

text(xoffset + nn + nn2 / 2, yoffset + mm / 2, ...
sprintf('%d', thisrk), ...
'HorizontalAlignment', 'center');
Expand All @@ -61,6 +87,20 @@ function spy_draw_block(H, xoffset, yoffset, rk)
thisrk = rank(H.B21);

fill(xb, yb, offdiag_color);

if strcmp(variant, 'svd')
% Draw the SVD plot of the block in the background
s = log(svd(H.B21));
nsvd = length(s);
for j = 1 : nsvd
sjm = (s(1) - s(j)) / (s(1) - s(end)) * mm2;
xbj = [xoffset + nn * (j-1)/nsvd, xoffset + nn * j / nsvd, ...
xoffset + nn * j / nsvd, xoffset + nn * (j-1)/nsvd ];
ybj = [yoffset + mm + sjm, yoffset + mm + sjm, yoffset + mm + mm2, yoffset + mm + mm2 ];
fill(xbj, ybj, svd_color, 'LineStyle','none');
end
end

text(xoffset + nn / 2, yoffset + mm + mm2 / 2, ...
sprintf('%d', thisrk), ...
'HorizontalAlignment', 'center');
Expand Down

0 comments on commit c0e4dbe

Please sign in to comment.