Remove Blocks of Matrix from a Matrix

Hello Guys,

I am actually working on a way to remove a certain block of a matrix A with given indices from the matrix A.

The problem is the following:

I have a symmetric matrix A which the size is (N*axis, N*axis), and I have a set of Nodes to Remove (NodesRemove) and a set of Nodes to keep (NodesKeep).

However, what I am seeking to perform is not to remove the line and column of the index, but to remove its complete block contribution (composed of axis line and axis column), simply put, from (idx_node*axis)-(axis-1) to (idx_node*axis) line and column.

To solve this, I have created a function "ReducedSystemParam" that can perform this operation perfectly, this works well and fast enough when N is small, however it is slow when N is bigger. (e.g. N >= 400)

The complete code is as follows:

clear;
close;
clc;

N            =  400;
axis         =  3;

A            =  rand(N*axis);
A            =  A - tril(A,-1) + triu(A,1)';

all_N        =  [1:N];
NodesRemove  =  [12; 32; 68; 95; 121; 145; 175; 182; 206; 242; 265; 287; 310; 335; 358; 389]; 
NodesKeep    =  (setdiff(all_N, NodesRemove))';

tic
[ Ar ]       =  ReducedSystemParam ( A, NodesKeep, axis );
toc

function [ Kr] = ReducedSystemParam ( K, NodesFree, Nbr_Axis )

Size_Kr             =  ( size ( NodesFree, 1 ) ) * Nbr_Axis;
Kr                  =  zeros(Size_Kr, Size_Kr);

% Extract Kr from K
for i = 1 : size (NodesFree, 1)
    for j = i : size (NodesFree, 1)      
        start_line_index_K     =  ( NodesFree(i) * Nbr_Axis ) - (Nbr_Axis - 1);
        finish_line_index_K    =  ( NodesFree(i) * Nbr_Axis );   
        start_column_index_K   =  ( NodesFree(j) * Nbr_Axis ) - (Nbr_Axis - 1);
        finish_column_index_K  =  ( NodesFree(j) * Nbr_Axis );
 
        start_line_index_Kr     =  ( i * Nbr_Axis ) - (Nbr_Axis - 1);
        finish_line_index_Kr    =  ( i * Nbr_Axis );   
        start_column_index_Kr   =  ( j * Nbr_Axis ) - (Nbr_Axis - 1);
        finish_column_index_Kr  =  ( j * Nbr_Axis );
        
        Kr ( start_line_index_Kr : finish_line_index_Kr , start_column_index_Kr : finish_column_index_Kr )  =  K ( start_line_index_K : finish_line_index_K , start_column_index_K : finish_column_index_K );
        Kr ( start_column_index_Kr : finish_column_index_Kr , start_line_index_Kr : finish_line_index_Kr )  =  ( Kr ( start_line_index_Kr : finish_line_index_Kr , start_column_index_Kr : finish_column_index_Kr ) )'; 
        
    end
end

end

Does anyone have any suggestions of how we can better solve this problem.

Thanks in advance

Answer 1

[Vidhi Agarwal]

Hi Amehri,

The following approach demonstrates an optimized version of the function:

By Incorporating ‘vectorization’ and ‘direct submatrix extraction’ significantly enhances the code's efficiency. Here is the explanation:

Vectorization Enhancement: The updated function utilizes vectorized operations through ‘arrayfun’ with 'UniformOutput' set to false. This approach efficiently transforms ‘NodesFree’ into exact row and column indices in one step, eliminating the need for nested loops and improving code efficiency.

Optimization via Direct Submatrix Extraction: This involves the direct utilization of pre-calculated indices, referred to as ‘actualIndicesKept’. With these indices, the function extracts the desired submatrix in one cohesive operation. This method not only simplifies the process but also leverages MATLAB's optimized matrix indexing and submatrix extraction functionalities to their fullest extent.

Following is the modified code using these strategies:

clear; 
close all; 
clc; 
N = 400; 
axis = 3; 
A = rand(N*axis); 
A = A - tril(A,-1) + triu(A,1)'; 
all_N = 1:N; 
NodesRemove = [12; 32; 68; 95; 121; 145; 175; 182; 206; 242; 265; 287; 310; 335; 358; 389];  
NodesKeep = setdiff(all_N, NodesRemove)'; 
tic 
Ar = ReducedSystemParamOptimized(A, NodesKeep, axis); 
toc 
function Kr = ReducedSystemParamOptimized(K, NodesFree, Nbr_Axis) 
    % Pre-calculate the indices based on NodesFree and Nbr_Axis 
    indices = arrayfun(@(x) (x-1)*Nbr_Axis+1 : x*Nbr_Axis, NodesFree, 'UniformOutput', false); 
    indices = horzcat(indices{:}); % Convert cell array to a vector of indices 
    % Directly extract the submatrix based on calculated indices 
    Kr = K(indices, indices); 
end 

For better understanding of ‘vectorization’ you can refer to this document:

• https://www.mathworks.com/help/matlab/matlab_prog/vectorization.html