正如here所暗示的，使用图形方法应该是可能的。下面是链接答案中描述的过程的基本实现：

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

让我们看看这对问题中的问题是如何起作用的：

% 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
%}

旁注：如果所有矩阵的大小相同，可以使用pagesvd函数（R2021b中引入）方便地将此计算矢量化，此函数不需要CPU计算工具箱。