gamma = 0.9
dim = 4

M = np.zeros((dim,dim))
for i in range(dim)[::-1]:
    M += np.diagflat([gamma**(dim-i-1)]*(i+1),dim-i-1)

print(M)

array([[ 1.   ,  0.9  ,  0.81 ,  0.729],
       [ 0.   ,  1.   ,  0.9  ,  0.81 ],
       [ 0.   ,  0.   ,  1.   ,  0.9  ],
       [ 0.   ,  0.   ,  0.   ,  1.   ]])