# Understand numpy indexing


```python
#%matplotlib widget
%matplotlib inline
%load_ext autoreload
%autoreload 2
```


```python
import numpy as np
import matplotlib.pyplot as plt
import sympy
```

## A few ways to get test numpy arrays


```python
np.arange(3), np.arange(4,8), np.arange(5,1,-2)
```




    (array([0, 1, 2]), array([4, 5, 6, 7]), array([5, 3]))



For experiments with multiplication, arrays of primes may be helpful:


```python
def arangep(n, starting_index=0):
    sympy.sieve.extend_to_no(starting_index + n)
    return np.array(sympy.sieve._list[starting_index:starting_index + n])
```


```python
arangep(5), arangep(4,2)
```




    (array([ 2,  3,  5,  7, 11]), array([ 5,  7, 11, 13]))



# Shapes and Indexing

Indexing [basics](https://numpy.org/devdocs/user/basics.indexing.html#basics-indexing) and [details](https://numpy.org/devdocs/reference/arrays.indexing.html#arrays-indexing)


```python
a = np.arange(2*3*4).reshape(2,3,4); print(a)
```

    [[[ 0  1  2  3]
      [ 4  5  6  7]
      [ 8  9 10 11]]
    
     [[12 13 14 15]
      [16 17 18 19]
      [20 21 22 23]]]


Indexing is row-major order (smallest-address-delta last) (C-style):


```python
a[0,0,1], a[0,1,0], a[1,0,0]
```




    (1, 4, 12)




```python
a[0], a[0,0], a[0,0,0]
```




    (array([[ 0,  1,  2,  3],
            [ 4,  5,  6,  7],
            [ 8,  9, 10, 11]]),
     array([0, 1, 2, 3]),
     0)




```python
a[0], a[0][0], a[0][0][0]
```




    (array([[ 0,  1,  2,  3],
            [ 4,  5,  6,  7],
            [ 8,  9, 10, 11]]),
     array([0, 1, 2, 3]),
     0)




```python
a.flat[7:12]
```




    array([ 7,  8,  9, 10, 11])



# Multiplicative-type operations


```python
a = arangep(2)
b = arangep(2,2)
a,b
```




    (array([2, 3]), array([5, 7]))



Binary scalar operations on vectors just map


```python
a+1, a*2, a+b, a*b, b/a, b%a
```




    (array([3, 4]),
     array([4, 6]),
     array([ 7, 10]),
     array([10, 21]),
     array([2.5       , 2.33333333]),
     array([1, 1]))



[`dot`](https://numpy.org/devdocs/reference/generated/numpy.dot.html) is "alternative matrix product with different broadcasting rules"


```python
a.dot(b), b.dot(a)
```




    (31, 31)




```python
m = arangep(4,4).reshape(2,2); m
```




    array([[11, 13],
           [17, 19]])



## Dot product

Matrix dot vector produces vector of dot products of rows of the matrix with the vector:


```python
m.dot(a), a.dot(m[0]), a.dot(m[1]), m[0], m[1]
```




    (array([61, 91]), 61, 91, array([11, 13]), array([17, 19]))



vector dot matrix produces vector of dot products of columns of matrix with the vector:


```python
a.dot(m), a.dot(m[:,0]), a.dot(m[:,1]), m[:,0], m[:,1]
```




    (array([73, 83]), 73, 83, array([11, 17]), array([13, 19]))



`@` is infix [matrix multiplication](https://numpy.org/devdocs/reference/generated/numpy.matmul.html#numpy.matmul)


```python
a, m, m @ a, a @ m, m.T @ a
```




    (array([2, 3]),
     array([[11, 13],
            [17, 19]]),
     array([61, 91]),
     array([73, 83]),
     array([73, 83]))



Right-multiplication by a matrix is equivalent to left-multiplication by its transpose:


```python
a @ m, m.T @ a, a @ m.T, m @ a
```




    (array([73, 83]), array([73, 83]), array([61, 91]), array([61, 91]))



### "Vectorizing" the dot product
e.g. when we batch inputs to the network. \
Imagine `a` and `b` are both to be run through a network which does multiplication by `m`


```python
c = 2*a + b
a, b, c, a @ m, b @ m, c @ m
```




    (array([2, 3]),
     array([5, 7]),
     array([ 9, 13]),
     array([73, 83]),
     array([174, 198]),
     array([320, 364]))



The convenient representation *(see below)*, is for the input vectors to be contiguous and adjacent in memory, as would happen if you read them into a memoryview of an array, and reshaped it appropriately, e.g.:


```python
X = np.array([2,3, 5,7, 9,13]).reshape(-1, 2); X
```




    array([[ 2,  3],
           [ 5,  7],
           [ 9, 13]])




```python
X @ m
```




    array([[ 73,  83],
           [174, 198],
           [320, 364]])




```python
X @ m + np.array([1000, 2000])
```




    array([[1073, 2083],
           [1174, 2198],
           [1320, 2364]])




```python
X.shape
```




    (3, 2)



## Einstein summation notation

Numpy provides [Einstein summation](https://mathworld.wolfram.com/EinsteinSummation.html) operations with [einsum](https://numpy.org/devdocs/reference/generated/numpy.einsum.html)
1. Repeated indices are implicitly summed over.
1. Each index can appear at most twice in any term.
1. Each term must contain identical non-repeated indices.


```python
es = np.einsum
```

 $$a_{ik}a_{ij} \equiv \sum_{i} a_{ik}a_{ij}$$

$$M_{ij}v_j=\sum_{j}M_{ij}v_j$$


```python
es('ij,j', m, a), es('ij,i', m, a)
```




    (array([61, 91]), array([73, 83]))




```python
es('j,ij', a, m), es('i,ij', a, m)
```




    (array([61, 91]), array([73, 83]))



Scalar multiplication bei


```python
all(es('ij,j', m, a) == es('j,ij', a, m))
```




    True



### Lorem Ipsum


```python
m2 = np.zeros((2,3), np.int); m2
```




    array([[0, 0, 0],
           [0, 0, 0]])




```python
m2[1] = np.arange(3); m2
```




    array([[0, 0, 0],
           [0, 1, 2]])




```python
m3 = arangep(8).reshape(4,2).T; m3
```




    array([[ 2,  5, 11, 17],
           [ 3,  7, 13, 19]])




```python
m3[:,0]
```




    array([2, 3])




```python
m @ m3[:,0]
```




    array([61, 91])




```python
h = m @ m3; h
```




    array([[ 61, 146, 290, 434],
           [ 91, 218, 434, 650]])




```python
b, b[...,np.newaxis]
```




    (array([5, 7]),
     array([[5],
            [7]]))




```python
h + b[...,np.newaxis]
```




    array([[ 66, 151, 295, 439],
           [ 98, 225, 441, 657]])



## Convenient representations

Suppose you have many __x__ to run through a net. What is the convenient representation?

Consider a two-input net, e.g. the XOR net. We want to vectorize the evaluation of the net, and its backprop. In the case of XOR the entire input domain is four vectors: { (0,0), (0,1), (1,0), (1,1) }:


```python
X = np.array([0,0, 0,1, 1,0, 1,1]).reshape(-1,2); X
```




    array([[0, 0],
           [0, 1],
           [1, 0],
           [1, 1]])



This is a convenient ordering for input, with each input vector contiguous in memory. But it's not in the form of column vectors for the classical left-multiplication by a transformation matrix to yield a column matrix product.


```python
m = np.arange(4).reshape(2,2) + 1; m
```




    array([[1, 2],
           [3, 4]])




```python
m @ np.array([1, 2]).reshape(2,1)
```




    array([[ 5],
           [11]])



We can transpose the input before left-multiplying ...


```python
m @ X.T
```




    array([[0, 2, 1, 3],
           [0, 4, 3, 7]])



... and transpose it back:


```python
Y = (m @ X.T).T; Y
```




    array([[0, 0],
           [2, 4],
           [1, 3],
           [3, 7]])



Or we can be less pedantic about expressing the matrix multiply:


```python
X @ m.T
```




    array([[0, 0],
           [2, 4],
           [1, 3],
           [3, 7]])



In Einstein summation notation:


```python
es('ij,kj', X, m)
```




    array([[0, 0],
           [2, 4],
           [1, 3],
           [3, 7]])



If we really require the matrix on the left, we can index thus:


```python
es('ij,kj->ki', m, X)
```




    array([[0, 0],
           [2, 4],
           [1, 3],
           [3, 7]])



---
### What way is faster?


```python
timeit(X @ m.T)
```

    1.05 µs ± 5.03 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)



```python
timeit(es('ij,kj->ki', m, X))
```

    2.16 µs ± 9.84 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)



```python
tm = m.T
```


```python
timeit(X @ tm)
```

    897 ns ± 9.11 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)


No surprise, fastest is to have the transposed matrix ready. No surprise that the Einstein summation is slower, as it requires formulating a loop from the string of indexes. But what if the input data is much larger? E.g.


```python
Xlarge = np.arange(2*10000).reshape(10000,2); Xlarge
```




    array([[    0,     1],
           [    2,     3],
           [    4,     5],
           ...,
           [19994, 19995],
           [19996, 19997],
           [19998, 19999]])




```python
timeit(Xlarge @ tm)
```

    86.9 µs ± 142 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)



```python
timeit(es('ij,kj->ki', m, Xlarge))
```

    126 µs ± 126 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)


The parsing of the index string and formulating a plan is maybe 1.6 µs, but the loop is 


```python
(156 + 1.4 - 3.01)/94.2
```




    1.63895966029724



64% slower.

---

Adding another vector to each result vector of the multiply:


```python
a, a + Y, Y + a
```




    (array([2, 3]),
     array([[ 2,  3],
            [ 4,  7],
            [ 3,  6],
            [ 5, 10]]),
     array([[ 2,  3],
            [ 4,  7],
            [ 3,  6],
            [ 5, 10]]))



Applying a function to each result:


```python
relu = np.vectorize(lambda x: max(0,x))
```

Try it out:


```python
t = arangep(10).reshape(5,2) - 12; t
```




    array([[-10,  -9],
           [ -7,  -5],
           [ -1,   1],
           [  5,   7],
           [ 11,  17]])




```python
relu(t)
```




    array([[ 0,  0],
           [ 0,  0],
           [ 0,  1],
           [ 5,  7],
           [11, 17]])



---


```python
X @ m.T
```




    array([[0, 0],
           [2, 4],
           [1, 3],
           [3, 7]])




```python
es('ij,kj', X, m)
```




    array([[0, 0],
           [2, 4],
           [1, 3],
           [3, 7]])




```python
X, m
```




    (array([[0, 0],
            [0, 1],
            [1, 0],
            [1, 1]]),
     array([[1, 2],
            [3, 4]]))



___

### Outer product


```python
a, b = arangep(2), arangep(3,2)
a, b
```




    (array([2, 3]), array([ 5,  7, 11]))




```python
es('i,j', a, b), es('j,i', a, b), np.outer(a, b), np.outer(b, a)
```




    (array([[10, 14, 22],
            [15, 21, 33]]),
     array([[10, 15],
            [14, 21],
            [22, 33]]),
     array([[10, 14, 22],
            [15, 21, 33]]),
     array([[10, 15],
            [14, 21],
            [22, 33]]))




```python
a, b = arangep(4).reshape(2,2), arangep(4,4).reshape(2,2)
a,b
```




    (array([[2, 3],
            [5, 7]]),
     array([[11, 13],
            [17, 19]]))



$$ \sum_j{outer(a[:,j],b[:,j])}$$


```python
es('...i,j', a, b)
```

# Vectorized dot product


```python
t = np.arange(4*3).reshape(-1,3)
t
```




    array([[ 0,  1,  2],
           [ 3,  4,  5],
           [ 6,  7,  8],
           [ 9, 10, 11]])




```python
[t[i].dot(t[i]) for i in range(4)]
```




    [5, 50, 149, 302]




```python
es('...i,...i', t, t)
```




    array([  5,  50, 149, 302])




```python
es('...i,...i', t[1:], t[:-1])
```




    array([ 14,  86, 212])



`np.arccos` range is $[0,\pi)$. Want to convert to $[-\pi/2, \pi/2)$.


```python
a = np.arange(-5,6)*np.pi/5
ca = np.cos(a)
aca = np.arccos(ca)
a, ca, aca
```




    (array([-3.14159265, -2.51327412, -1.88495559, -1.25663706, -0.62831853,
             0.        ,  0.62831853,  1.25663706,  1.88495559,  2.51327412,
             3.14159265]),
     array([-1.        , -0.80901699, -0.30901699,  0.30901699,  0.80901699,
             1.        ,  0.80901699,  0.30901699, -0.30901699, -0.80901699,
            -1.        ]),
     array([3.14159265, 2.51327412, 1.88495559, 1.25663706, 0.62831853,
            0.        , 0.62831853, 1.25663706, 1.88495559, 2.51327412,
            3.14159265]))




```python
b = a >= (np.pi/2)
b
```




    array([False, False, False, False, False,  True,  True,  True,  True,
            True])




```python
b * np.pi/2
```




    array([0.        , 0.        , 0.        , 0.        , 0.        ,
           1.57079633, 1.57079633, 1.57079633, 1.57079633, 1.57079633])




```python
a[a>=np.pi/2] -= np.pi
a
```




    array([ 0.        ,  0.31415927,  0.62831853,  0.9424778 ,  1.25663706,
           -1.57079633, -1.25663706, -0.9424778 , -0.62831853, -0.31415927])




```python
a = np.arange(8)
a, a%3, a%3>0
```




    (array([0, 1, 2, 3, 4, 5, 6, 7]),
     array([0, 1, 2, 0, 1, 2, 0, 1]),
     array([False,  True,  True, False,  True,  True, False,  True]))




```python
a%3>0
```

# END
---


```python
t =  np.array([0,0, 0,1, 1,0, 1,1]).reshape(4,2); t
```


```python
f = lambda a, b: 1 if (a > 0.5) ^ (b > 0.5) else 0
```


```python
f(1,0), f(1,1)
```


```python
f(t[1,0], t[1,1])
```


```python
[f(x[0], x[1]) for x in t]
```


```python
def exor(a, b):
    return 1 if (a > 0.5) ^ (b > 0.5) else 1
```


```python
#np.vectorize(exor, signature='(i)->()')(t)
```


```python
f2 = lambda v: 1 if (v[0] > 0.5) ^ (v[1] > 0.5) else 0
```


```python
#np.vectorize(f2)(t)
```


```python
np.vectorize(f2, signature='(i)->()')(t)
```


```python
[a for a in t[0]]
```

---


```python
a = np.arange(25).reshape(5,5)
b = np.arange(5)
c = np.arange(6).reshape(2,3)
```


```python
a,b,c
```


```python
np.einsum('ii', a)
```


```python
np.einsum('ii->i', a)
```


```python
np.trace(a)
```


```python
np.einsum('ji', a)
```


```python
np.einsum('ji,i', a, b)
```


```python
a.dot(b)
```


```python
a[:,0]
```


```python
a[:,0].dot(b)
```


```python
a[:,1]
```


```python
d = np.arange(125).reshape(5,5,5)
```


```python
np.einsum('iii', d)
```


```python
sum([d[i][i][i] for i in range(5)])
```


```python
np.einsum('iij',d)
```


```python
np.einsum('iiz', d)
```


```python
[sum([d[i][i][j] for i in range(5)]) for j in range(5)]
```


```python
sum(a[:])
```


```python
a[0]
```


```python
timeit(np.einsum('iii', d))
```


```python
es = np.einsum
```


```python
es('ijk,kji',d,d)
```


```python
timeit(es('iii', d))
```


```python
es('i,ij', b, a)
```


```python
es('ij', a)
```


```python
es('i', b)
```


```python
g = np.arange(4).reshape(2,2)
```


```python
g
```


```python
es('ij,jk',g,g)
```


```python
g@g
```


```python
g[:]
```


```python
h = np.arange(2); h
```


```python
h.dot(g)
```


```python
es('i,ij', h, g)
```


```python
g.dot(h)
```


```python
es('ji,i', g, h)
```


```python
es('ij,j', g, h)
```


```python
g[0,1]
```


```python
np.array(1)
```


```python
np.array([1])
```


```python
np.array(2)
```


```python
np.array([1,2])
```


```python
np.array([2])
```


```python
np.array(0)
```


```python
np.array(0).shape
```


```python
np.array([0]).shape
```


```python
np.array(0)+1
```


```python
np.array([0])+1
```


```python
np.array(3).dot(np.array(5))
```


```python
np.array([3]).dot(np.array(5))
```


```python
np.array([3]).dot(np.array([5]))
```


```python
es('i,i', np.array([3]), np.array([5]))
```


```python
#es('i,i', np.array(3), np.array(5))
```


```python
es('', np.array(3), np.array(5))
```


```python
int(np.array(3))
```


```python
int(np.array([3]))
```


```python
a = np.arange(4).reshape(2,2) + 1; print(a)
```


```python
b = np.arange(2) + 1; print(b)
```


```python
b.dot(a)
```


```python
es('ij, jk', a, np.array([[1,0],[1,0]]))
```


```python
a[:,0]
```


```python
a[:,1]
```


```python
sum(a[:,0])
```


```python
a.T
```


```python
es('...j->...', a)
```


```python
a,b
```


```python
a.dot(b), b.dot(a)
```


```python
b.shape
```


```python
c = b.reshape(2,1); c
```


```python
b.dot(c), b@c
```


```python
es('...i,i...', b, c)
```


```python
timeit(b.dot(c))
```


```python
timeit(b@c)
```


```python
timeit(es('...i,i...', b, c))
```


```python
timeit(es('i,i...', b, c))
```


```python
a,b,c
```


```python
a@b
```


```python
a@c
```


```python
b.shape, c.shape
```


```python
es('ij,j...', a,c)
```


```python
es('ij,j', a, b)
```


```python
es('ij,j...', a, b)
```


```python
es('ij,j->i', a, b)
```


```python
es('ij,j->j', a, b)
```


```python
es('ij,i...', a,c)
```


```python
es('ij,j...', a, a)
```


```python
es('...j,ij', a, a)
```


```python
Xd = np.array([0,0,1,0,0,1,1,1]).reshape(2,4); Xd
```


```python
a
```


```python
a@Xd
```


```python
Xd.reshape(4,2)
```


```python
t = np.arange(8).reshape(4,2); t
```


```python
a @ t.T
```


```python
t.T
```


```python
Ellipsis
```


```python
b
```


```python
b[:, np.newaxis]
```


```python
t
```


```python
t[np.newaxis]
```


```python
x = np.arange(3); x
```


```python
x[:,np.newaxis] + x[np.newaxis,:]
```


```python
x[:,np.newaxis] * x[np.newaxis,:]
```


```python
x[np.newaxis,:] * x[:,np.newaxis]
```


```python
x[:,np.newaxis], x[np.newaxis,:]
```


```python
t.dot(np.arange(2) + 1)
```


```python
(np.arange(2) + 1).dot(t)
```


```python
t,a
```


```python
t.dot(a)
```


```python
np.prime
```


```python
import sympy
```


```python
np.array(list(sympy.sieve.primerange(2000,2050)))
```


```python
pa = lambda n: np.array([sympy.prime(i+1) for i in range(n)])
```


```python
pa(50)
```


```python
tp = pa(1000)
```


```python
tp[-1]
```


```python
tp2 = pa(10000)
```


```python
np.sign(np.arange(5)-2)
```


```python
np.array([1]*2).shape
```


```python
np.ones(np.array([3,5]).shape)
```


```python
M=np.arange(4).reshape(2,2)
b=np.arange(2)+1
x=np.arange(2)+5
```


```python
M@x + b, (M@x) + b
```


```python
a, b = arangep(2), arangep(3,2)
a, b
```




    (array([2, 3]), array([ 5,  7, 11]))




```python
t = es('i,j', a, b)
t
```




    array([[10, 14, 22],
           [15, 21, 33]])




```python
t[0].dot(t[0]), t[1].dot(t[1]), sum([t[0].dot(t[0]), t[1].dot(t[1])]), es('ij,ij', t, t)
```




    (780, 1755, 2535, 2535)




```python
#es('ij,ij', t, t)
np.atleast_2d(t)
```




    array([[10, 14, 22],
           [15, 21, 33]])




```python

```

85a67d7a3e645990736550a1aac47c1acca4d991
ipynb
Jupyter Notebook
nbs/matrix_operations.ipynb
pramasoul/aix
98333b875f6c6cda6dee86e6eab02c5ddc622543
MIT
nbs/matrix_operations.ipynb
pramasoul/aix
98333b875f6c6cda6dee86e6eab02c5ddc622543
MIT
2021-11-29T03:44:00.000Z
2021-12-19T05:34:04.000Z
nbs/matrix_operations.ipynb
pramasoul/aix
98333b875f6c6cda6dee86e6eab02c5ddc622543
MIT
Qwen/Qwen-72B
1. YES
2. YES