Why does calling np.array() on this list comprehension produce a 3d array instead of 2d?

Question

I have a script produces the first several iterations of a Markov matrix multiplying a given set of input values. With the matrix stored as A and the start values in the column u0, I use this list comprehension to store the output in an array:

out = np.array([ ( (A**n) * u0).T for n in range(10) ])

The output has shape (10,1,6), but I want the output in shape (10,6) instead. Obviously, I can fix this with .reshape(), but is there a way to avoid creating the extra dimension in the first place, perhaps by simplifying the list comprehension or the inputs?

Here's the full script and output:

import numpy as np

# Random 6x6 Markov matrix
n = 6
A = np.matrix([ (lambda x: x/x.sum())(np.random.rand(n)) for _ in range(n)]).T
print(A)

#[[0.27457312 0.20195133 0.14400801 0.00814027 0.06026188 0.23540134]
# [0.21526648 0.17900277 0.35145882 0.30817386 0.15703758 0.21069114]
# [0.02100412 0.05916883 0.18309142 0.02149681 0.22214047 0.15257011]
# [0.17032696 0.11144443 0.01364982 0.31337906 0.25752732 0.1037133 ]
# [0.03081507 0.2343255  0.2902935  0.02720764 0.00895182 0.21920371]
# [0.28801424 0.21410713 0.01749843 0.32160236 0.29408092 0.07842041]]

# Random start values
u0 = np.matrix(np.random.randint(51, size=n)).T
print(u0)

#[[31]
# [49]
# [44]
# [29]
# [10]
# [ 0]]

# Find the first 10 iterations of the Markov process
out = np.array([ ( (A**n) * u0).T for n in range(10) ])
print(out)

#[[[31.         49.         44.         29.         10.
#    0.        ]]
#
# [[25.58242101 41.41600236 14.45123543 23.00477134 26.08867045
#   32.45689942]]
#
# [[26.86917065 36.02438292 16.87560159 26.46418685 22.66236879
#   34.10428921]]
#
# [[26.69224394 37.06346073 16.59208202 26.48817955 22.56696872
#   33.59706504]]
#
# [[26.68772374 36.99727159 16.49987315 26.5003184  22.61130862
#   33.7035045 ]]
#
# [[26.68766363 36.98517264 16.50532933 26.51717543 22.592951
#   33.71170797]]
#
# [[26.68695152 36.98895204 16.50314718 26.51729716 22.59379049
#   33.70986161]]
#
# [[26.68682195 36.98848867 16.50286371 26.51763013 22.59362679
#   33.71056876]]
#
# [[26.68681128 36.98850409 16.50286036 26.51768807 22.59359453
#   33.71054167]]
#
# [[26.68680313 36.98851046 16.50285038 26.51769497 22.59359219
#   33.71054886]]]

print(out.shape)

#(10, 1, 6)

out = out.reshape(10,n)
print(out)

#[[31.         49.         44.         29.         10.          0.        ]
# [25.58242101 41.41600236 14.45123543 23.00477134 26.08867045 32.45689942]
# [26.86917065 36.02438292 16.87560159 26.46418685 22.66236879 34.10428921]
# [26.69224394 37.06346073 16.59208202 26.48817955 22.56696872 33.59706504]
# [26.68772374 36.99727159 16.49987315 26.5003184  22.61130862 33.7035045 ]
# [26.68766363 36.98517264 16.50532933 26.51717543 22.592951   33.71170797]
# [26.68695152 36.98895204 16.50314718 26.51729716 22.59379049 33.70986161]
# [26.68682195 36.98848867 16.50286371 26.51763013 22.59362679 33.71056876]
# [26.68681128 36.98850409 16.50286036 26.51768807 22.59359453 33.71054167]
# [26.68680313 36.98851046 16.50285038 26.51769497 22.59359219 33.71054886]]

Answer 1

I think your confusion lies with how arrays can be joined.

Start with a simple 1d array (in numpy 1d is a real thing, not just a 'row vector' or 'column vector'):

In [288]: arr = np.arange(6)                                                    
In [289]: arr                                                                   
Out[289]: array([0, 1, 2, 3, 4, 5])

np.array joins element arrays along a new 1st dimension:

In [290]: np.array([arr,arr])                                                   
Out[290]: 
array([[0, 1, 2, 3, 4, 5],
       [0, 1, 2, 3, 4, 5]])

np.stack with the default axis value does the same thing. Read its docs.

We can make a 2d array, a column vector:

In [291]: arr1 = arr[:,None]                                                    
In [292]: arr1                                                                  
Out[292]: 
array([[0],
       [1],
       [2],
       [3],
       [4],
       [5]])
In [293]: arr1.shape                                                            
Out[293]: (6, 1)

Using np.array on its transpose the (1,6) arrays:

In [294]: np.array([arr1.T, arr1.T])                                            
Out[294]: 
array([[[0, 1, 2, 3, 4, 5]],

       [[0, 1, 2, 3, 4, 5]]])
In [295]: _.shape                                                               
Out[295]: (2, 1, 6)

Note the middle size 1 dimension, that bothered you.

np.vstack joins the arrays along the existing 1st dimension. It does not add one:

In [296]: np.vstack([arr1.T, arr1.T])                                           
Out[296]: 
array([[0, 1, 2, 3, 4, 5],
       [0, 1, 2, 3, 4, 5]])

Or we could join the arrays horizontally, on the 2nd dimension:

In [297]: np.hstack([arr1, arr1])                                               
Out[297]: 
array([[0, 0],
       [1, 1],
       [2, 2],
       [3, 3],
       [4, 4],
       [5, 5]])

That is (6,2) which can be transposed to (2,6):

In [298]: np.hstack([arr1, arr1]).T                                             
Out[298]: 
array([[0, 1, 2, 3, 4, 5],
       [0, 1, 2, 3, 4, 5]])

Answer 2

I made a few changes to the code, although I'm not 100% certain that the result is still the same (I am not familiar with Markov chains).

import numpy as np

n = 6

num_proc_iters = 10

rand_nums_arr = np.random.random_sample((n, n))
rand_nums_arr = np.transpose(rand_nums_arr / rand_nums_arr.sum(axis=1))

u0 = np.random.randint(51, size=n)

res_arr = np.concatenate([np.linalg.matrix_power(rand_nums_arr, curr) @ u0 for curr in range(num_proc_iters)])

I would love to hear if anyone can think of any further improvements.

Answer 3

If you use np.array() for input and @ for matrix multiplication, it works as expected.

# Random 6x6 Markov matrix
n = 6
A = np.array([ (lambda x: x/x.sum())(np.random.rand(n)) for _ in range(n)]).T

# Random start values
u0 = np.random.randint(51, size=n).T

# Find the first 10 iterations of the Markov process
out = np.array([ ( np.linalg.matrix_power(A,n) @ u0).T for n in range(10) ])
print(out)

#[[29.         24.          5.         12.         10.         32.        ]
# [15.82875119 13.53436868 20.61648725 19.22478172 20.34082205 22.45478912]
# [21.82434718 10.06037119 14.29281935 20.75271393 18.76134538 26.30840297]
# [20.77484848 10.1379821  15.47488423 19.4965479  20.05618311 26.05955418]
# [21.02944236 10.09401438 15.24263478 19.48662616 19.95767996 26.18960236]
# [20.96887722 10.11647819 15.30729334 19.44261102 20.00089222 26.16384802]
# [20.98086362 10.11522779 15.29529799 19.44899285 19.99137187 26.16824587]
# [20.97795615 10.11606978 15.29817734 19.44798612 19.99293494 26.16687566]
# [20.97858032 10.11591954 15.29752865 19.44839852 19.99245389 26.16711909]
# [20.97844343 10.11594666 15.29766432 19.4483417  19.99254284 26.16706104]]