I want to perform SVD on M separate problems of potentially different sizes. The sizes of the datasets (matrices) to decompose are N_M⨯2, where N_M is different for each of the M problems and ∀M, N_M > 2.
An obvious way to do this is through a {potentially parallelized} loop, where each dataset is fed into SVD separately and the results are handled separately as well.
Another possibility is to arrange the different datasets as diagonal blocks in a {potentially sparse} matrix, then call SVD exactly once, thereby vectorizing the computation and having to deal with only a single set of outputs. This is the approach I would like to focus on in the question.
The common implementations of SVD (e.g. MATLAB's and Numpy's) sort the outputs (columns/rows) of the U and V matrices according to the magnitude of singular values (S). This auto-sorting is detrimental for our block-diagonal case since the association between the singular values/vectors and the blocks that generated them is lost. For example,
%% Define matrices to decompose:
A = [1, 1, 2;
2, 1, 3];
B = [1, 2;
2, 1];
M = blkdiag(A,B);
%{
M =
1 1 2 0 0
2 1 3 0 0
0 0 0 1 2
0 0 0 2 1
%}
%% Apply SVD
[Ua, Sa, Va] = svd(A, 'econ');
[Ub, Sb, Vb] = svd(B, 'econ');
[Um, Sm, Vm] = svd(M, 'econ');
%{
Ua = # Ub = # Um =
-0.5449 -0.8385 # -0.7071 -0.7071 # -0.5449 0 0 -0.8385
-0.8385 0.5449 # -0.7071 0.7071 # -0.8385 0 0 0.5449
# # 0 -0.7071 0.7071 0
# # 0 -0.7071 -0.7071 0
# #
Sa = # Sb = # Sm =
4.4552 0 # 3.0000 0 # 4.4552 0 0 0
0 0.3888 # 0 1.0000 # 0 3.0000 0 0
# # 0 0 1.0000 0
# # 0 0 0 0.3888
# #
Va = # Vb = # Vm =
-0.4987 0.6465 # -0.7071 0.7071 # -0.4987 0 0 0.6465
-0.3105 -0.7551 # -0.7071 -0.7071 # -0.3105 0 0 -0.7551
-0.8092 -0.1087 # # -0.8092 0 0 -0.1087
# # 0 -0.7071 -0.7071 0
# # 0 -0.7071 0.7071 0
%}
Note how Um, Sm and Vm have the same values as {Ua, Ub}, {Sa, Sb} and {Va, Vb}, accordingly, but the values originating from A are not grouped together in any of the decompositions of M.
If so, how can we rearrange ("unsort") Um, Sm, and Vm such that they correspond to the input blocks? After reorganization, the expected result is:
% If M = blkdiag(A, B), we want (up to sign indeterminacy):
Um == blkdiag(Ua, Ub);
Sm == blkdiag(Sa, Sb);
Vm == blkdiag(Va, Vb);
As hinted here, this should be possible using a graph approach. Below is a basic implementation of the procedure described in the linked answer:
function [row_perms, col_perms, PL, PR] = block_diagonalize(A)
% Find the permutation of rows/columns of any rectangular array A such
% that it becomes block-diagonal.
%
% Assumptions:
% 1) The input is block-diagonalizeable.
% 2) The input has no all-zero rows or columns.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Find connected components
nnz = sparse(A ~= 0);
[nr, nc] = size(A);
connectivity = [zeros(nr, 'logical'), nnz;
nnz.', zeros(nc, 'logical')];
G = graph(connectivity, ["r" + (1:nr), "c" + (1:nc)]);
[comp, ~] = G.conncomp; % compare: comp = conncomp(G,'OutputForm','cell');
%% TODO: Compare the 2nd output of `conncomp` to a new "block size" input
%% TODO: Add handling of zero rows/columns to the ordering
if max(comp) > min(nr, nc)
warning("No all-zero rows/columns allowed! An additional " + ...
"invocation of this function will be needed.")
end
%% Get the sorting permutations:
row_order = comp(1:nr);
col_order = comp(nr+1:end);
if issorted(row_order) && issorted(col_order)
% No reordering necessary
row_perms = row_order;
col_perms = col_order;
return
end
[~, row_perms] = sort(row_order);
[~, col_perms] = sort(col_order);
%% Build the bidirectional permutation matrices
if nargout > 2
PL = sparse(1:nr, row_perms, ones(nr, 1), nr, nr);
PR = sparse(col_perms, 1:nc, ones(nc, 1), nc, nc);
% PL * A * PR
end
end
Let's see how this works for the problem in the question:
% We assume that [Um, Sm, Vm] have already been computed.
[urp, ucp, upl, upr] = block_diagonalize(Um);
[vrp, vcp, vpl, vpr] = block_diagonalize(Vm);
%{
>> upl * Um * upr
ans =
-0.5449 -0.8385 0 0
-0.8385 0.5449 0 0
0 0 -0.7071 0.7071
0 0 -0.7071 -0.7071
>> vpl * Vm * vpr
ans =
-0.4987 0.6465 0 0
-0.3105 -0.7551 0 0
-0.8092 -0.1087 0 0
0 0 -0.7071 -0.7071
0 0 -0.7071 0.7071
>> assert(isequal(ucp, vcp)); d = diag(Sm); diag(d(ucp))
ans =
4.4552 0 0 0
0 0.3888 0 0
0 0 3.0000 0
0 0 0 1.0000
>> norm(M - (upl * Um * upr)*diag(d(ucp))*(vpl * Vc * vpr).')
ans =
8.8818e-16
%}
Sidenote: if the sizes of all matrices are the same, the pagesvd function (introduced in R2021b) can be used to conveniently vectorize this computation. This function requires no toolboxes for CPU computations.