Create 4D upper diagonal array from 3D
<p>Let's say that I have a <code>(x, y, z)</code> sized matrix. Now, I wish to create a new matrix of dimension <code>(x, y, i, i)</code>, where the <code>(i, i)</code> matrix is upper diagonal and constructed from the values on the <code>z</code>-dimension. Is there some easy way of doing this in <code>numpy</code> without using more than 1 for-loop (looping over x)? Thanks.</p>

<p><strong>EDIT</strong></p>

<pre><code>original = np.array([
    [
        [0, 1, 3],
        [4, 5, 6]
    ],
    [
        [7, 8, 9],
        [3, 2, 1]
    ],
])

new = np.array([
    [
        [
            [0, 1],
            [0, 3]
        ],
        [
            [4, 5],
            [0, 6]
        ]
    ],
    [
        [
            [7, 8],
            [0, 9]
        ],
        [
            [3, 2],
            [0, 1]
        ]
    ]
])
</code></pre>

<p>So, using the above we see that</p>

<pre><code>original[0, 0, :] = [0 1 3]

new[0, 0, :, :] = [[0 1]
                  [0 3]]
</code></pre>

<p>Here's an approach using <code>boolean-indexing</code> -</p>

<pre><code>n = 2 # This would depend on a.shape[-1]
out = np.zeros(a.shape[:2] + (n,n,),dtype=a.dtype)
out[:,:,np.arange(n)[:,None] &lt;= np.arange(n)] = a
</code></pre>

<p>Sample run -</p>

<pre><code>In [247]: a
Out[247]: 
array([[[0, 1, 3],
        [4, 5, 6]],

       [[7, 8, 9],
        [3, 2, 1]]])

In [248]: out
Out[248]: 
array([[[[0, 1],
         [0, 3]],

        [[4, 5],
         [0, 6]]],


       [[[7, 8],
         [0, 9]],

        [[3, 2],
         [0, 1]]]])
</code></pre>

<p>Another approach could be suggested using <code>subscripted-indexing</code> to replace the last step -</p>

<pre><code>r,c = np.triu_indices(n)
out[:,:,r,c] = a
</code></pre>

<p><strong>Note :</strong> As stated earlier, <code>n</code> would depend on <code>a.shape[-1]</code>. Here, we had <code>a.shape[-1]</code> as <code>3</code>, so <code>n</code> was <code>2</code>. If <code>a.shape[-1]</code> were <code>6</code>, <code>n</code> would be <code>3</code> and so on. The relationship is : <code>(n*(n+1))//2 == a.shape[-1]</code>.</p>

python numpy vectorization