Skip to article content

Pre-pre-school

Back to Article
Chapter 2: Vectors, matrices and multidimensional arrays
Download Notebook

Chapter 2: Vectors, matrices and multidimensional arrays

import numpy as np

The NumPy array object

data = np.array([[1, 2], [3, 4], [5, 6]])
type(data)
numpy.ndarray
data
array([[1, 2], [3, 4], [5, 6]])
data.ndim
2
data.shape
(3, 2)
data.size
6
data.dtype
dtype('int64')
data.nbytes
48

Data types

np.array([1, 2, 3], dtype=int)
array([1, 2, 3])
np.array([1, 2, 3], dtype=float)
array([1., 2., 3.])
np.array([1, 2, 3], dtype=complex)
array([1.+0.j, 2.+0.j, 3.+0.j])
data = np.array([1, 2, 3], dtype=float)
data
array([1., 2., 3.])
data.dtype
dtype('float64')
data = np.array([1, 2, 3], dtype=int)
data.dtype
dtype('int64')
data
array([1, 2, 3])
data = np.array([1, 2, 3], dtype=float)
data
array([1., 2., 3.])
data.astype(int)
array([1, 2, 3])
d1 = np.array([1, 2, 3], dtype=float)
d2 = np.array([1, 2, 3], dtype=complex)
d1 + d2
array([2.+0.j, 4.+0.j, 6.+0.j])
(d1 + d2).dtype
dtype('complex128')
np.sqrt(np.array([-1, 0, 1]))
/tmp/ipykernel_25889/208196152.py:1: RuntimeWarning: invalid value encountered in sqrt
  np.sqrt(np.array([-1, 0, 1]))
array([nan, 0., 1.])
np.sqrt(np.array([-1, 0, 1], dtype=complex))
array([0.+1.j, 0.+0.j, 1.+0.j])

Real and imaginary parts

data = np.array([1, 2, 3], dtype=complex)
data
array([1.+0.j, 2.+0.j, 3.+0.j])
data.real
array([1., 2., 3.])
data.imag
array([0., 0., 0.])

Creating arrays

Arrays created from lists and other array-like objects

np.array([1, 2, 3, 4])
array([1, 2, 3, 4])
data.ndim
1
data.shape
(3,)
np.array([[1, 2], [3, 4]])
array([[1, 2], [3, 4]])
data.ndim
1
data.shape
(3,)

Arrays filled with constant values

np.zeros((2, 3))
array([[0., 0., 0.], [0., 0., 0.]])
np.ones(4)
array([1., 1., 1., 1.])
data = np.ones(4)
data.dtype
dtype('float64')
data = np.ones(4, dtype=np.int64)
data.dtype
dtype('int64')
x1 = 5.4 * np.ones(10)
x2 = np.full(10, 5.4)
x1 = np.empty(5)
x1.fill(3.0)
x1
array([3., 3., 3., 3., 3.])
x2 = np.full(5, 3.0)
x2
array([3., 3., 3., 3., 3.])

Arrays filled with incremental sequences

np.arange(0.0, 11, 1)
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])
np.linspace(0, 10, 11)
array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.])

Arrays filled with logarithmic sequences

np.logspace(0, 2, 5)  # 5 data points between 10**0=1 to 10**2=100
array([ 1. , 3.16227766, 10. , 31.6227766 , 100. ])

Mesh-grid arrays

x = np.array([-1, 0, 1])
y = np.array([-2, 0, 2])
X, Y = np.meshgrid(x, y)
X
array([[-1, 0, 1], [-1, 0, 1], [-1, 0, 1]])
Y
array([[-2, -2, -2], [ 0, 0, 0], [ 2, 2, 2]])
Z = (X + Y) ** 2
Z
array([[9, 4, 1], [1, 0, 1], [1, 4, 9]])

Creating uninitialized arrays

np.empty(3, dtype=float)
array([1., 2., 3.])

Creating arrays with properties of other arrays

def f(x):
    y = np.ones_like(x)
    # compute with x and y
    return y

Creating matrix arrays

np.identity(4)
array([[1., 0., 0., 0.], [0., 1., 0., 0.], [0., 0., 1., 0.], [0., 0., 0., 1.]])
np.eye(3, k=1)
array([[0., 1., 0.], [0., 0., 1.], [0., 0., 0.]])
np.eye(3, k=-1)
array([[0., 0., 0.], [1., 0., 0.], [0., 1., 0.]])
np.diag(np.arange(0, 20, 5))
array([[ 0, 0, 0, 0], [ 0, 5, 0, 0], [ 0, 0, 10, 0], [ 0, 0, 0, 15]])

Index and slicing

One-dimensional arrays

a = np.arange(0, 11)
a
array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10])
a[0]  # the first element
np.int64(0)
a[-1]  # the last element
np.int64(10)
a[4]  # the fifth element, at index 4
np.int64(4)
a[1:-1]
array([1, 2, 3, 4, 5, 6, 7, 8, 9])
a[1:-1:2]
array([1, 3, 5, 7, 9])
a[:5]
array([0, 1, 2, 3, 4])
a[-5:]
array([ 6, 7, 8, 9, 10])
a[::-2]
array([10, 8, 6, 4, 2, 0])

Multidimensional arrays

f = lambda m, n: n + 10 * m
A = np.fromfunction(f, (6, 6), dtype=int)
A
array([[ 0, 1, 2, 3, 4, 5], [10, 11, 12, 13, 14, 15], [20, 21, 22, 23, 24, 25], [30, 31, 32, 33, 34, 35], [40, 41, 42, 43, 44, 45], [50, 51, 52, 53, 54, 55]])
A[:, 1]  # the second column
array([ 1, 11, 21, 31, 41, 51])
A[1, :]  # the second row
array([10, 11, 12, 13, 14, 15])
A[:3, :3]  # upper half diagonal block matrix
array([[ 0, 1, 2], [10, 11, 12], [20, 21, 22]])
A[3:, :3]  # lower left off-diagonal block matrix
array([[30, 31, 32], [40, 41, 42], [50, 51, 52]])
A[::2, ::2]  # every second element starting from 0, 0
array([[ 0, 2, 4], [20, 22, 24], [40, 42, 44]])
A[1::2, 1::3]  # every second element starting from 1, 1
array([[11, 14], [31, 34], [51, 54]])

Views

B = A[1:5, 1:5]
B
array([[11, 12, 13, 14], [21, 22, 23, 24], [31, 32, 33, 34], [41, 42, 43, 44]])
B[:, :] = 0
A
array([[ 0, 1, 2, 3, 4, 5], [10, 0, 0, 0, 0, 15], [20, 0, 0, 0, 0, 25], [30, 0, 0, 0, 0, 35], [40, 0, 0, 0, 0, 45], [50, 51, 52, 53, 54, 55]])
C = B[1:3, 1:3].copy()
C
array([[0, 0], [0, 0]])
C[:, :] = 1  # this does not affect B since C is a copy of the view B[1:3, 1:3]
C
array([[1, 1], [1, 1]])
B
array([[0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]])

Fancy indexing and Boolean-valued indexing

A = np.linspace(0, 1, 11)
A[np.array([0, 2, 4])]
array([0. , 0.2, 0.4])
A[[0, 2, 4]]
array([0. , 0.2, 0.4])
A > 0.5
array([False, False, False, False, False, False, True, True, True, True, True])
A[A > 0.5]
array([0.6, 0.7, 0.8, 0.9, 1. ])
A = np.arange(10)
A
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
indices = [2, 4, 6]
B = A[indices]
B[0] = -1  # this does not affect A
A
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
A[indices] = -1
A
array([ 0, 1, -1, 3, -1, 5, -1, 7, 8, 9])
A = np.arange(10)
B = A[A > 5]
B[0] = -1  # this does not affect A
A
array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
A[A > 5] = -1
A
array([ 0, 1, 2, 3, 4, 5, -1, -1, -1, -1])

Reshaping and resizing

data = np.array([[1, 2], [3, 4]])
np.reshape(data, (1, 4))
array([[1, 2, 3, 4]])
data.reshape(4)
array([1, 2, 3, 4])
data = np.array([[1, 2], [3, 4]])
data
array([[1, 2], [3, 4]])
data.flatten()
array([1, 2, 3, 4])
data.flatten().shape
(4,)
data = np.arange(0, 5)
column = data[:, np.newaxis]
column
array([[0], [1], [2], [3], [4]])
row = data[np.newaxis, :]
row
array([[0, 1, 2, 3, 4]])
data = np.arange(5)
data
array([0, 1, 2, 3, 4])
np.vstack((data, data, data))
array([[0, 1, 2, 3, 4], [0, 1, 2, 3, 4], [0, 1, 2, 3, 4]])
data = np.arange(5)
data
array([0, 1, 2, 3, 4])
np.hstack((data, data, data))
array([0, 1, 2, 3, 4, 0, 1, 2, 3, 4, 0, 1, 2, 3, 4])
data = data[:, np.newaxis]
np.hstack((data, data, data))
array([[0, 0, 0], [1, 1, 1], [2, 2, 2], [3, 3, 3], [4, 4, 4]])

Vectorized expressions

Arithmetic operations

x = np.array([[1, 2], [3, 4]])
y = np.array([[5, 6], [7, 8]])
x + y
array([[ 6, 8], [10, 12]])
y - x
array([[4, 4], [4, 4]])
x * y
array([[ 5, 12], [21, 32]])
y / x
array([[5. , 3. ], [2.33333333, 2. ]])
x * 2
array([[2, 4], [6, 8]])
2**x
array([[ 2, 4], [ 8, 16]])
y / 2
array([[2.5, 3. ], [3.5, 4. ]])
(y / 2).dtype
dtype('float64')
x = np.array([1, 2, 3, 4]).reshape(2, 2)
z = np.array([1, 2, 3, 4])
# x / z
z = np.array([[2, 4]])
z
array([[2, 4]])
z.shape
(1, 2)
x / z
array([[0.5, 0.5], [1.5, 1. ]])
zz = np.concatenate([z, z], axis=0)
zz
array([[2, 4], [2, 4]])
x / zz
array([[0.5, 0.5], [1.5, 1. ]])
z = np.array([[2], [4]])
z.shape
(2, 1)
x / z
array([[0.5 , 1. ], [0.75, 1. ]])
zz = np.concatenate([z, z], axis=1)
zz
array([[2, 2], [4, 4]])
x / zz
array([[0.5 , 1. ], [0.75, 1. ]])
x = np.array([[1, 3], [2, 4]])
x = x + y
x
array([[ 6, 9], [ 9, 12]])
x = np.array([[1, 3], [2, 4]])
x += y
x
array([[ 6, 9], [ 9, 12]])

Elementwise functions

x = np.linspace(-1, 1, 11)
x
array([-1. , -0.8, -0.6, -0.4, -0.2, 0. , 0.2, 0.4, 0.6, 0.8, 1. ])
y = np.sin(np.pi * x)
np.round(y, decimals=4)
array([-0. , -0.5878, -0.9511, -0.9511, -0.5878, 0. , 0.5878, 0.9511, 0.9511, 0.5878, 0. ])
np.add(np.sin(x) ** 2, np.cos(x) ** 2)
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
np.sin(x) ** 2 + np.cos(x) ** 2
array([1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.])
def heaviside(x):
    return 1 if x > 0 else 0
heaviside(-1)
0
heaviside(1.5)
1
x = np.linspace(-5, 5, 11)
# heaviside(x)
heaviside = np.vectorize(heaviside)
heaviside(x)
array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
def heaviside(x):
    return 1.0 * (x > 0)

Aggregate functions

data = np.random.normal(size=(15, 15))
np.mean(data)
np.float64(0.09261224693203207)
data.mean()
np.float64(0.09261224693203207)
data = np.random.normal(size=(5, 10, 15))
data.sum(axis=0).shape
(10, 15)
data.sum(axis=(0, 2)).shape
(10,)
data.sum()
np.float64(6.048928142352301)
data = np.arange(1, 10).reshape(3, 3)
data
array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
data.sum()
np.int64(45)
data.sum(axis=0)
array([12, 15, 18])
data.sum(axis=1)
array([ 6, 15, 24])

Boolean arrays and conditional expressions

a = np.array([1, 2, 3, 4])
b = np.array([4, 3, 2, 1])
a < b
array([ True, True, False, False])
np.all(a < b)
np.False_
np.any(a < b)
np.True_
if np.all(a < b):
    print("All elements in a are smaller than their corresponding element in b")
elif np.any(a < b):
    print("Some elements in a are smaller than their corresponding elemment in b")
else:
    print("All elements in b are smaller than their corresponding element in a")
Some elements in a are smaller than their corresponding elemment in b

x = np.array([-2, -1, 0, 1, 2])
x > 0
array([False, False, False, True, True])
1 * (x > 0)
array([0, 0, 0, 1, 1])
x * (x > 0)
array([0, 0, 0, 1, 2])
def pulse(x, position, height, width):
    return height * (x >= position) * (x <= (position + width))
x = np.linspace(-5, 5, 11)
pulse(x, position=-2, height=1, width=5)
array([0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0])
pulse(x, position=1, height=1, width=5)
array([0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1])
def pulse(x, position, height, width):
    return height * np.logical_and(x >= position, x <= (position + width))
x = np.linspace(-4, 4, 9)
np.where(x < 0, x**2, x**3)
array([16., 9., 4., 1., 0., 1., 8., 27., 64.])
np.select([x < -1, x < 2, x >= 2], [x**2, x**3, x**4])
array([ 16., 9., 4., -1., 0., 1., 16., 81., 256.])
np.choose([0, 0, 0, 1, 1, 1, 2, 2, 2], [x**2, x**3, x**4])
array([ 16., 9., 4., -1., 0., 1., 16., 81., 256.])
x[abs(x) > 2]
array([-4., -3., 3., 4.])
np.nonzero(abs(x) > 2)
(array([0, 1, 7, 8]),)
x[np.nonzero(abs(x) > 2)]
array([-4., -3., 3., 4.])

Set operations

a = np.unique([1, 2, 3, 3])
b = np.unique([2, 3, 4, 4, 5, 6, 5])
np.isin(a, b)
array([False, True, True])
1 in a
True
1 in b
False
np.all(np.isin(a, b))
np.False_
np.union1d(a, b)
array([1, 2, 3, 4, 5, 6])
np.intersect1d(a, b)
array([2, 3])
np.setdiff1d(a, b)
array([1])
np.setdiff1d(b, a)
array([4, 5, 6])

Operations on arrays

data = np.arange(9).reshape(3, 3)
data
array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
np.transpose(data)
array([[0, 3, 6], [1, 4, 7], [2, 5, 8]])
data = np.random.randn(1, 2, 3, 4, 5)
data.shape
(1, 2, 3, 4, 5)
data.T.shape
(5, 4, 3, 2, 1)

Matrix and vector operations

A = np.arange(1, 7).reshape(2, 3)
A
array([[1, 2, 3], [4, 5, 6]])
B = np.arange(1, 7).reshape(3, 2)
B
array([[1, 2], [3, 4], [5, 6]])
np.dot(A, B)
array([[22, 28], [49, 64]])
A @ B
array([[22, 28], [49, 64]])
np.dot(B, A)
array([[ 9, 12, 15], [19, 26, 33], [29, 40, 51]])
B @ A
array([[ 9, 12, 15], [19, 26, 33], [29, 40, 51]])
A = np.arange(9).reshape(3, 3)
A
array([[0, 1, 2], [3, 4, 5], [6, 7, 8]])
x = np.arange(3)
x
array([0, 1, 2])
np.dot(A, x)
array([ 5, 14, 23])
A.dot(x)
array([ 5, 14, 23])
A = np.random.rand(3, 3)
B = np.random.rand(3, 3)
Ap = B @ A @ np.linalg.inv(B)
Ap
array([[ 0.95824093, -5.59107895, 3.66628953], [ 1.12937933, -6.01346248, 3.41382176], [ 2.14182688, -11.37408349, 6.35465289]])
Ap = np.dot(B, np.dot(A, np.linalg.inv(B)))
Ap
array([[ 0.95824093, -5.59107895, 3.66628953], [ 1.12937933, -6.01346248, 3.41382176], [ 2.14182688, -11.37408349, 6.35465289]])
Ap = B.dot(A.dot(np.linalg.inv(B)))
Ap
array([[ 0.95824093, -5.59107895, 3.66628953], [ 1.12937933, -6.01346248, 3.41382176], [ 2.14182688, -11.37408349, 6.35465289]])
A = np.matrix(A)
B = np.matrix(B)
Ap = B * A * B.I
A = np.asmatrix(A)
B = np.asmatrix(B)
Ap = B * A * B.I
Ap = np.asarray(Ap)
Ap
array([[ 0.95824093, -5.59107895, 3.66628953], [ 1.12937933, -6.01346248, 3.41382176], [ 2.14182688, -11.37408349, 6.35465289]])
np.inner(x, x)
np.int64(5)
np.dot(x, x)
np.int64(5)
y = x[:, np.newaxis]
y
array([[0], [1], [2]])
np.dot(y.T, y)
array([[5]])
x = np.array([1, 2, 3])
np.outer(x, x)
array([[1, 2, 3], [2, 4, 6], [3, 6, 9]])
np.kron(x, x)
array([1, 2, 3, 2, 4, 6, 3, 6, 9])
np.kron(x[:, np.newaxis], x[np.newaxis, :])
array([[1, 2, 3], [2, 4, 6], [3, 6, 9]])
np.kron(np.ones((2, 2)), np.identity(2))
array([[1., 0., 1., 0.], [0., 1., 0., 1.], [1., 0., 1., 0.], [0., 1., 0., 1.]])
np.kron(np.identity(2), np.ones((2, 2)))
array([[1., 1., 0., 0.], [1., 1., 0., 0.], [0., 0., 1., 1.], [0., 0., 1., 1.]])
x = np.array([1, 2, 3, 4])
y = np.array([5, 6, 7, 8])
np.einsum("n,n", x, y)
np.int64(70)
np.inner(x, y)
np.int64(70)
A = np.arange(9).reshape(3, 3)
B = A.T
np.einsum("mk,kn", A, B)
array([[ 5, 14, 23], [ 14, 50, 86], [ 23, 86, 149]])
np.all(np.einsum("mk,kn", A, B) == np.dot(A, B))
np.True_
References
  1. Johansson, R. (2024). Numerical Python: Scientific Computing and Data Science Applications with Numpy, SciPy and Matplotlib. Apress. 10.1007/979-8-8688-0413-7
Numericalpython
Chapter 1: Computing with Python
Numericalpython
Chapter 2: Figures