the dimension of matrix when using 'sparse' to generate sparse matrix

Hi@Dyuman Joshi, Thank you for your feedback. My answer adds clarification on why users might think 3D worked, a code demonstration, and a concrete solution (ndSparse) for true N-dimensional sparse arrays. These details provide unique value. In the spirit of community, it’s always more helpful to contribute clarifying thoughts rather than just flagging. So,could you please share your thoughts to help @Li over here, I would be more happy to hear what you have to say.

@Li - MATLAB's tril() and sparse() functions have never supported 3D arrays, so nothing was actually removed between 2017b and 2025a. I would like a second opinion from @Walter and @Torsten who are experts and do respect their opinions always in the post. It might be you had the ndSparse toolbox installed back then and forgot about it. For your specific problem, you need to apply tril to each 8×8 page individually using a loop, and for the block diagonal creation, the issue with H(8,8,:) is that it extracts only single elements (giving you 1×1×20) rather than the full pages—you need to use num2cell(H, [1 2]) to convert each page to a cell, then pass that to blkdiag. I've put together working code (took some time to implement to help you out) in the artifact below that includes tril3d() and blkdiag3d() utility functions that handle this correctly, and everything works fine with complex doubles (the built-in blkdiag has always supported complex matrices, so if you're getting an error about that, check "which blkdiag -all" to see if a custom function is shadowing the built-in). The code demonstrates the correct approach for your 8×8×20 case and creates the sparse block diagonal matrix you need.

Here is the code

% Solution for Li's 3D matrix problem
% Problem: tril() doesn't work on 3D arrays, block diagonal syntax 
%issues
%% Your immediate problem - applying tril to 3D matrices

% your data (8x8 pages, 20 deep, complex)
H = rand(8,8,20) + 1i*rand(8,8,20);
fprintf('Starting with %dx%dx%d complex matrix\n', size(H,1), 
size(H,2),   
size(H,3));

% tril doesn't accept 3D - you have to loop through pages
% there's no built-in function for this unfortunately
H_lower = zeros(size(H), 'like', H);  % pre-allocate (keeps it 
complex)
for i = 1:size(H,3)
  H_lower(:,:,i) = tril(H(:,:,i));
end

fprintf('Applied tril to all %d pages\n', size(H,3));

% quick check that it worked
if nnz(triu(H_lower(:,:,1), 1)) == 0
  fprintf('Looks good - upper triangles are zero\n');
end

%% Block diagonal - THIS is probably where your syntax was wrong

% WRONG WAY (what you might have tried):
% BD = blkdiag(H(8,8,:))  % <- this extracts single elements!

% The problem: H(8,8,:) gives you a 1x1x20 array of just the (8,8)  
%elements You need the full 8x8 pages

% RIGHT WAY:
pages = num2cell(H_lower, [1 2]);  % converts each page to a cell
BD = blkdiag(pages{:});             % now this works correctly

fprintf('\nBlock diagonal created: %dx%d\n', size(BD,1), size(BD,2));

% convert to sparse if you want
BD_sparse = sparse(BD);
sparsity = 100 * (1 - nnz(BD_sparse)/numel(BD_sparse));
fprintf('Sparse version is %.1f%% sparse\n', sparsity);

% verify it's still complex
if ~isreal(BD_sparse)
  fprintf('Complex values preserved: %.3f%+.3fi (sample)\n', ...
      real(full(BD_sparse(1,1))), imag(full(BD_sparse(1,1))));
end

%% Utility functions you can reuse

% I made these into functions so you don't have to write loops every 
%time

function result = tril3d(A, k)
  % apply tril to each page of a 3D array
  % usage: result = tril3d(A) or result = tril3d(A, k) for diagonal 
  %offset

    if nargin < 2
        k = 0;  % main diagonal by default
    end

    result = zeros(size(A), 'like', A);
    for i = 1:size(A,3)
        result(:,:,i) = tril(A(:,:,i), k);
    end
  end

function result = triu3d(A, k)
  % same thing but for upper triangular

    if nargin < 2
        k = 0;
    end

    result = zeros(size(A), 'like', A);
    for i = 1:size(A,3)
        result(:,:,i) = triu(A(:,:,i), k);
    end
  end

function BD = blkdiag3d(A, make_sparse)
  % creates block diagonal from 3D array pages
  % usage: BD = blkdiag3d(A) or BD = blkdiag3d(A, true) for sparse

    if nargin < 2
        make_sparse = false;
    end

    pages = num2cell(A, [1 2]);
    BD = blkdiag(pages{:});

    if make_sparse
        BD = sparse(BD);
    end
  end

%% Test the utility functions

fprintf('\n--- Testing utility functions ---\n');

test_data = rand(4,4,5) + 1i*rand(4,4,5);

% test tril3d
lower = tril3d(test_data);
fprintf('tril3d: processed %dx%dx%d -> check upper triangle: ',   
size(test_data));
if nnz(triu(lower(:,:,1),1)) == 0
  fprintf('pass\n');
else
  fprintf('fail\n');
end

% test triu3d  
upper = triu3d(test_data);
fprintf('triu3d: processed %dx%dx%d -> check lower triangle: ',   
size(test_data));
if nnz(tril(upper(:,:,1),-1)) == 0
  fprintf('pass\n');
else
  fprintf('fail\n');
end

% test blkdiag3d
bd_test = blkdiag3d(test_data);
fprintf('blkdiag3d: created %dx%d block diagonal\n', size(bd_test));

% with sparse option
bd_sparse_test = blkdiag3d(test_data, true);
fprintf('blkdiag3d (sparse): %d nonzeros out of %d elements\n', ...
  nnz(bd_sparse_test), numel(bd_sparse_test));

%% Complete workflow example with your dimensions

fprintf('\n--- Complete example with your 8x8x20 case ---\n');

% generate data
H_example = rand(8,8,20) + 1i*rand(8,8,20);

% apply lower triangular
H_example_lower = tril3d(H_example);

% create sparse block diagonal
BD_example = blkdiag3d(H_example_lower, true);

fprintf('Input: %dx%dx%d complex matrix\n', size(H_example));
fprintf('Output: %dx%d sparse matrix\n', size(BD_example));
fprintf('Nonzeros: %d (%.2f%% full)\n', nnz(BD_example), ...
  100*nnz(BD_example)/numel(BD_example));

% check memory usage
mem_full = 8 * 160 * 160 * 2;  % bytes for complex double
mem_sparse = whos('BD_example').bytes;
fprintf('Memory: full would be %.1f MB, sparse is %.1f MB\n', ...
  mem_full/1e6, mem_sparse/1e6);

%% Common mistakes to avoid

fprintf('\n--- Common mistakes ---\n');

H_mistake = rand(8,8,20);

% Mistake 1: wrong indexing
wrong = H_mistake(8,8,:);
fprintf('H(8,8,:) gives size %dx%dx%d <- extracts SINGLE ELEMENTS\n',   
size(wrong));

correct = num2cell(H_mistake, [1 2]);
fprintf('num2cell(H,[1 2]) gives %d cells of %dx%d each <- THIS is   
what you want\n', ...
length(correct), size(correct{1}));

% Mistake 2: trying tril on 3D directly
fprintf('\nDirect tril(3D_array) will fail:\n');
try
  bad = tril(H_mistake);
  fprintf('  No error? Unexpected.\n');
catch ME
  fprintf('  Error: %s\n', ME.message);
  fprintf('  Solution: use tril3d() or loop through pages\n');
end

% Mistake 3: sparse on 3D
fprintf('\nDirect sparse(3D_array) will fail:\n');
try
  bad = sparse(H_mistake);
  fprintf('  No error? Unexpected.\n');
catch ME
  fprintf('  Error: %s\n', ME.message);
  fprintf('  Solution: create block diagonal first, then sparse 
it\n');
end

fprintf('\n--- Done ---\n');
fprintf('Copy the utility functions (tril3d, triu3d, blkdiag3d) to use
in your code\n');

Results: please see attached

Let me know how it goes.

Hi @Li,Now I see exactly what's happening - your error screenshots confirm you're using the ndSparse toolbox, and specifically its tril() function is throwing the error because it has a hard-coded check that refuses to work on arrays with more than 2 dimensions (even though ndSparse is designed for N-dimensional arrays, which seems like a bug or limitation in that version). This explains why it worked before - you likely had a different version of ndSparse or your code was structured slightly differently. Here are two solutions you can use right away:

Option 1 - Quick fix with your current ndSparse setup (just replace your line)

% Your original line that's failing:
% HL=spblkdiag(H1(:,:,1:end-1)+H1(:,:,2:end))+
 (-1i)*tril(spblkdiag(H2(:,:,1:end-1)+H2(:,:,2:end)));

% Replace with this (applies tril BEFORE spblkdiag):
H2_sum = H2(:,:,1:end-1) + H2(:,:,2:end);
H2_lower = zeros(size(H2_sum));
for k = 1:size(H2_sum, 3)
  H2_lower(:,:,k) = tril(H2_sum(:,:,k));
end
HL = spblkdiag(H1(:,:,1:end-1)+H1(:,:,2:end)) + 
(-1i)*spblkdiag(H2_lower);

Option 2 - Use pure MATLAB without ndSparse dependency (copy the utility functions from my previous post, then use)

H2_sum = H2(:,:,1:end-1) + H2(:,:,2:end);
H2_lower = tril3d(H2_sum);
HL = blkdiag3d(H1(:,:,1:end-1)+H1(:,:,2:end), true) + 
(-1i)*blkdiag3d(H2_lower, 
true);

Option 1 should work immediately with your current setup. Option 2 removes the ndSparse dependency entirely if you prefer that approach.

You are close solving your problem. You got this.

Let me know how it goes.

Hi @Li, Thankyou for your kind words. The reason your code worked before but doesn't now is likely one of these scenarios:

Most likely: You updated the ndSparse toolbox to a newer (or older) version that changed how tril() handles 3D arrays. The ndSparse tril() implementation has a dimension check that throws an error for arrays with more than 2 dimensions, even though ndSparse is designed for N-dimensional sparse arrays. This seems like an inconsistency in the toolbox itself.

Other possibilities that I can hypothetically think about:

1. Your H2 array structure was different before (maybe it was 2D and you reshaped/concatenated it differently)

2. You were using MATLAB's built-in tril() instead of ndSparse's tril() (path order changed)

3. Your previous code applied tril() before creating the 3D array

To confirm what changed, check these things:
% Check which tril you're using
which tril -all

% Check your ndSparse version (if it has version info)
ver ndSparse

% Check H2 dimensions
size(H2)
ndims(H2)

If you want to avoid the for-loop and keep your original one-liner style, you could modify ndSparse's tril.m file directly (if you have edit access) to handle 3D arrays, or create a wrapper function that handles the 3D case automatically.

Again, the for-loop solution I provided is actually quite efficient since it's just looping over the third dimension and MATLAB handles 2D tril() operations very quickly. For most applications, the performance difference will be negligible compared to the spblkdiag operations.

If you find your old code version and it works, compare which tril() function it's calling - that will tell you exactly what changed.