Python – Nested loops with a variable number of loops

I have a function that outputs the matrix on a new basis. However, depending on the size of the matrix, the number of base matrices varies. So in the simplified “Matlab pseudocode”:

if matrixsize==1
   for a1=1:4
      out(a1)=Matrix*basis(a1)
   end

elseif matrixsize==2
   for a1=1:4
      for a2=a1:4
         out(a1,a2)=Matrix*basis(a1)*basis(a2)
      end
   end

elseif matrixsize==3
   for a1=1:4
      for a2=a1:4
          for a3=a2:4
             out(a1,a2,a3)=Matrix*basis(a1)*basis(a2)*basis(a3)
          end
      end
   end

elseif ...

And so on

Is it possible to write this code for any value of the matrix size?
In other words: is it possible to create a loop that automatically creates the loop above?
If this doesn’t work in Matlab, is there a solution in Python?

(Background: This question comes from quantum physics, I want to write about the quantum state of Pauliki).

This is a valid Matlab code that shows the problem:

function T=newbasis(n)

%create a random matrix
m=2^n;
M=randn(m);

%Pauli matrices
s{1}=sparse([1,0;0,1]);
s{2}=sparse([0,1; 1,0]);
s{3}=sparse([0,-1i; 1i,0]);
s{4}=sparse([1,0;0,-1]);

if n==1
    for a1=1:4
        T(a1)=trace(M*betterkron(s{a1}));
    end

elseif n==2
    for a1=1:4
        for a2=a1:4
            T(a1,a2)=trace(M*betterkron(s{a1},s{a2}));
        end
    end

elseif n==3
    for a1=1:4
        for a2=a1:4
            for a3=a2:4    
                T(a1,a2,a3)=trace(M*betterkron(s{a1},s{a2},s{a3}));
            end
        end
    end    

else
    T=[]
end

%Not very clever but just to keep it simple
function krn=betterkron(A,varargin)
    krn = A;
    for j = 2:nargin;
        krn = kron(krn,varargin{j-1});
    end   
end

end

Although in principle such multiple loops can be done with recursive functions, it will be complicated. Fortunately, using multiple loops is not the best approach. MATLAB allows you to convert back and forth between N-dimensional subscripts and 1-dimensional linear indexes. So you can make a loop on the linear index and then convert it back to an N-dimensional subscript. So like this :

for i=1:numel(Matrix)  % loop over linear index 
    inds = ind2sub(size(Matrix), i);  % convert linear index to subscript

% Each index should be greater than or equal to the previous
    % e.g. a2=a1:4, a2 starts at a1 so cannot be less than a1
    if any(diff(inds) < 0)
        continue
    end

% Do the calculation
    % s{inds} is equivalent to s{i1}, s{i2}, ...
    T(i) = trace(M*betterkron(s{inds}));
end