Lesson 5

Overview

  • NumPy
  1. Initializing an array
  2. Initialization tricks - zeros/ones/empty
  • Specifying a datatype
  1. Shape and reshape
  2. Mathematical operations
  • Transposing multiple dimensions
  1. Accessing elements/indexing/slicing
  • Accessing elements/indexing
  • Slicing
  1. np.arange and np.linspace
  2. Tile and repeat
  3. Broadcasting
  • Cartesian products using broadcasting
  1. Random generators
  2. np.diff
  3. np.cumsum
  4. np.argmin/np.argmax/np.argsort
  5. np.where/np.nonzero
  6. np.bincount
  7. Dealing with big data
  • Sparse matrices
  • Matrix multiplication with diagonal matrices

NumPy

Vectors/arrays/matrices in Python

import numpy as np

1. Initializing an array

a = [
    [1, 2, 3],
    [4, 5, 6],
    [7, 8, 9]
]

Nested list representation

a
[[1, 2, 3], [4, 5, 6], [7, 8, 9]]

NumPy representation

b = np.array(a)
b
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])

2. Initialization tricks - zeros/ones/empty

c = np.zeros((5, 10))
c
array([[0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
       [0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]])
c = np.ones((5, 10))
c
array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])
c = np.empty((5, 10))
c
array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1.]])

2.1 Specifying a datatype

c = np.ones((5, 10), dtype=int)
c
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1]])

3. Shape and reshape

Access a matrix’s dimensions with .shape, and reshape its dimensions with .reshape()

c.shape
(5, 10)

Note: reshapes are NOT in-place

c.reshape((2, 25))
array([[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1],
       [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1]])
c.shape
(5, 10)

4. Mathematical operations

Let a denote a NumPy array, then: - + - element-wise addition - - - element-wise subtraction - * - element-wise multiplication - / - element-wise division - ** - element-wise powers - @ - matrix multiplication - np.linalg.inv(a) - matrix inversion - a.T or a.transpose() or np.transpose(a) - matrix transpose

Mathematical operations don’t work with lists

a * a
TypeError: can't multiply sequence by non-int of type 'list'

But they work with NumPy arrays

b * b
array([[ 1,  4,  9],
       [16, 25, 36],
       [49, 64, 81]])
b @ b
array([[ 30,  36,  42],
       [ 66,  81,  96],
       [102, 126, 150]])
np.linalg.inv(b) @ b
array([[ 14.,  12.,  10.],
       [-14.,  -8.,  -2.],
       [  0.,   0.,   0.]])

4.1 Transposing multiple dimensions

np.transpose(b, [1, 0])
array([[1, 4, 7],
       [2, 5, 8],
       [3, 6, 9]])

5. Accessing elements/indexing/slicing

5.1 Accessing elements/indexing

With lists, access each layer separately

a
[[1, 2, 3], [4, 5, 6], [7, 8, 9]]
print(a[0])
print(a[0][0])
[1, 2, 3]
1

With NumPy arrays this also works, but can also access multiple indices simultaneously

print(b[0])
print(b[0][0])
[1 2 3]
1
b[0, 0]
np.int64(1)

Can even access indices using tuples/lists/arrays

With a tuple, it is the same as if you typed in the index directly

idx = (0, 0)
b[idx]
np.int64(1)

With lists and arrays, the entire list indexes one dimension

idx = [0, 0]
b[idx]
array([[1, 2, 3],
       [1, 2, 3]])
b[[0, 0], :]
array([[1, 2, 3],
       [1, 2, 3]])
idx = np.array((0, 0))
b[idx]
array([[1, 2, 3],
       [1, 2, 3]])

Casting from list to tuple does work

idx = [0, 0]
b[tuple(idx)]
np.int64(1)

Accessing elements with tuples/lists/arrays doesn’t work for nested lists

a[idx]
TypeError: list indices must be integers or slices, not list

5.2 Slicing

Use a comma to separate dimensions and a colon (:) to access an entire dimension

a
[[1, 2, 3], [4, 5, 6], [7, 8, 9]]
a[:1]
[[1, 2, 3]]
b
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
b[:, [0, 2]]
array([[1, 3],
       [4, 6],
       [7, 9]])
b[:, [0, 2]][[0, 2], :]
array([[1, 3],
       [7, 9]])

You should probably avoid doing this, but you can choose to access only the first relevant dimensions of an array, and the remaining dimensions will implicitly have :

b[1, :]
array([4, 5, 6])
b[1]
array([4, 5, 6])

Note that slicing can allow you to update entire arrays in-place, saving allocation time or allowing in-place updates inside of functions

b[:, :] = b.copy()

6. np.arange and np.linspace

arange is the NumPy equivalent of Range, and linspace creates evenly-spaced values

c = np.arange(3)
c
array([0, 1, 2])
d = np.linspace(0, 1, 9)
d
array([0.   , 0.125, 0.25 , 0.375, 0.5  , 0.625, 0.75 , 0.875, 1.   ])

7. Tile and repeat

Tile stacks an entire array along a dimension (like tiles)

Repeat repeats each element of an array a given number of times, then repeats the next element, etc.

np.tile(c, 3)
array([0, 1, 2, 0, 1, 2, 0, 1, 2])
np.tile(c, (3, 2))
array([[0, 1, 2, 0, 1, 2],
       [0, 1, 2, 0, 1, 2],
       [0, 1, 2, 0, 1, 2]])
np.repeat(c, 3)
array([0, 0, 0, 1, 1, 1, 2, 2, 2])
np.repeat(c, (1, 2, 3))
array([0, 1, 1, 2, 2, 2])

8. Broadcasting

If certain dimensions between two arrays match, but one array has more dimensions than the other, NumPy can introduce ‘dummy’ dimensions into the smaller array so that they match

print(b.shape)
print(c.shape)
(3, 3)
(3,)
print(b)
print(c)
[[1 2 3]
 [4 5 6]
 [7 8 9]]
[0 1 2]
b + c
array([[ 1,  3,  5],
       [ 4,  6,  8],
       [ 7,  9, 11]])
(b.T + c).T
array([[ 1,  2,  3],
       [ 5,  6,  7],
       [ 9, 10, 11]])

Indexing with None creates a new dimension

print(c.shape)
print(c[None, :].shape)
(3,)
(1, 3)
b + c[None, :]
array([[ 1,  3,  5],
       [ 4,  6,  8],
       [ 7,  9, 11]])
b + c[:, None]
array([[ 1,  2,  3],
       [ 5,  6,  7],
       [ 9, 10, 11]])

8.1 Cartesian products using broadcasting

Use None in opposite dimensions

c
array([0, 1, 2])
c[:, None] + c[None, :]
array([[0, 1, 2],
       [1, 2, 3],
       [2, 3, 4]])

9. Random generators

Allows you to specify random seeds easily (but you don’t need to specify the seed)

seed = 1
rng = np.random.default_rng(seed)
rng.normal(size=10)
array([ 0.34558419,  0.82161814,  0.33043708, -1.30315723,  0.90535587,
        0.44637457, -0.53695324,  0.5811181 ,  0.3645724 ,  0.2941325 ])
rng.normal(size=10)
array([ 0.02842224,  0.54671299, -0.73645409, -0.16290995, -0.48211931,
        0.59884621,  0.03972211, -0.29245675, -0.78190846, -0.25719224])

Resetting the seed will replicate the results

seed = 1
rng = np.random.default_rng(seed)
rng.normal(size=10)
array([ 0.34558419,  0.82161814,  0.33043708, -1.30315723,  0.90535587,
        0.44637457, -0.53695324,  0.5811181 ,  0.3645724 ,  0.2941325 ])

10. np.diff

Take element-wise differences

a = rng.normal(size=100_000)
a
array([ 0.02842224,  0.54671299, -0.73645409, ..., -0.322686  ,
       -2.06155586, -0.70793319])
a[1:] - a[: -1]
array([ 0.51829075, -1.28316707,  0.57354414, ..., -0.84049398,
       -1.73886986,  1.35362267])
np.diff(a)
array([ 0.51829075, -1.28316707,  0.57354414, ..., -0.84049398,
       -1.73886986,  1.35362267])

11. np.cumsum

Cumulative sums

a = np.arange(5)
a
array([0, 1, 2, 3, 4])
np.cumsum(a)
array([ 0,  1,  3,  6, 10])

12. np.argmin/np.argmax/np.argsort

a = rng.normal(size=100)
a.argmax()
np.int64(79)

Note that argmax returns the first max

np.argmax([2, 2, 1])
np.int64(0)
a.argmin()
np.int64(3)
a.argsort()
array([ 3, 70, 66, 74, 31, 89, 65, 48, 68, 29, 30, 38, 96, 72,  7, 73, 43,
        2, 22, 33, 14,  0, 67, 34, 69, 95, 85, 83, 71, 28, 44, 77, 55, 15,
       19, 54, 49, 10, 64, 17, 36, 40,  9, 61, 46, 62, 47, 20, 88, 45, 21,
       81, 76,  5, 82, 60, 86, 58, 42, 25, 26, 97, 80, 53, 32, 11, 24, 41,
       13, 16, 87, 63, 12, 90, 98, 50, 78,  4, 37, 99,  8, 92, 56, 35, 52,
       18, 59, 84,  1, 39, 94, 51, 57, 23, 27, 75, 93, 91,  6, 79])
a[a.argsort()]
array([-2.87386515, -2.50281052, -2.44897667, -2.02623831, -1.96096236,
       -1.84160532, -1.74691692, -1.49858894, -1.48775121, -1.30809565,
       -1.2477091 , -1.22499879, -1.05225675, -1.01464596, -1.01122946,
       -0.99045975, -0.95871657, -0.91019895, -0.87470722, -0.84895448,
       -0.62521544, -0.61148934, -0.50171587, -0.49170749, -0.47036667,
       -0.45144132, -0.43364287, -0.41975742, -0.39152134, -0.39034982,
       -0.37303249, -0.36560617, -0.35555605, -0.34748432, -0.33698805,
       -0.32220437, -0.31917957, -0.20032613, -0.17710449, -0.16066618,
       -0.13862805, -0.12442159, -0.07698548, -0.04518314, -0.04500716,
       -0.02567932, -0.01250562, -0.00917651,  0.01075284,  0.01286459,
        0.02618942,  0.05272496,  0.06415227,  0.07496199,  0.09586831,
        0.10015544,  0.17405887,  0.18082354,  0.21079727,  0.24012108,
        0.28900234,  0.29361471,  0.38736811,  0.402049  ,  0.4096016 ,
        0.44171488,  0.44937894,  0.45799067,  0.50166376,  0.52507281,
        0.52898658,  0.54627952,  0.54725373,  0.63114246,  0.66714902,
        0.74472177,  0.77641005,  0.78511367,  0.84145412,  0.87439556,
        0.87609588,  0.91611766,  1.00685178,  1.03633602,  1.14481533,
        1.14522707,  1.19054368,  1.22460094,  1.25835347,  1.27815454,
        1.31836713,  1.37849169,  1.43487458,  1.46723572,  1.69429576,
        1.77365136,  1.80736638,  1.94029017,  1.95936278,  2.15130806])

13. np.where/np.nonzero

Briefly discussed in lecture 2 - find the indices where some condition holds

Supposedly np.nonzero is slightly faster

duplicate_max_lst = np.array([2, 2, 1])
np.argmax(duplicate_max_lst)
np.int64(0)
np.nonzero((duplicate_max_lst == duplicate_max_lst.max()).astype(int))[0]
array([0, 1])
np.where(duplicate_max_lst == duplicate_max_lst.max())[0]
array([0, 1])

14. np.bincount

Groupby-sum in NumPy. By default, it counts the number of observations in each group, but weights can be specified so it computes the sum of the weights within each group.

n = 100
n_ids = 10
ids = rng.integers(n_ids, size=n)
person_val = rng.uniform(size=n)
np.bincount(ids)
array([10,  8,  7, 13,  5, 13,  9, 15, 14,  6])

Can specify minlength to make sure the result of calling bincount is an array of the correct length - this helps in case the last group has no observations

np.bincount(ids, minlength=n_ids + 1)
array([10,  8,  7, 13,  5, 13,  9, 15, 14,  6,  0])
%%timeit
np.bincount(ids)
390 ns ± 14.4 ns per loop (mean ± std. dev. of 7 runs, 1,000,000 loops each)

bincount is much faster than Pandas groupby

import pandas as pd
id_df = pd.DataFrame({'id': ids, 'val': person_val})
id_df.groupby('id').size()
id
0    10
1     8
2     7
3    13
4     5
5    13
6     9
7    15
8    14
9     6
dtype: int64
%%timeit
id_df.groupby('id').size()
69.2 μs ± 1.29 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

bincount can use weights, so it is more like a groupby-sum

np.bincount(ids, person_val)
array([4.42350691, 2.9612293 , 4.11062489, 5.20176002, 2.56946445,
       7.53672527, 5.06821397, 6.52509301, 6.55220953, 2.85668455])
id_df.groupby('id')['val'].sum()
id
0    4.423507
1    2.961229
2    4.110625
3    5.201760
4    2.569464
5    7.536725
6    5.068214
7    6.525093
8    6.552210
9    2.856685
Name: val, dtype: float64

15. Dealing with big data

15.1 Sparse matrices

A sparse matrix is a matrix with mostly 0s. Since the only values that need to be known are the non-zero values, rather than storing the entire matrix, store only the non-zero entries.

from scipy.sparse import csc_matrix

The first argument specifies the values of the entries

The second argument is a tuple that specifies the index for each entry

The third argument specifies the shape

nvals = 5
J = csc_matrix((nvals - 1 - np.arange(nvals), (np.arange(nvals), nvals - 1 - np.arange(nvals))), shape=(nvals, nvals))
J.todense()
matrix([[0, 0, 0, 0, 4],
        [0, 0, 0, 3, 0],
        [0, 0, 2, 0, 0],
        [0, 1, 0, 0, 0],
        [0, 0, 0, 0, 0]])
J = csc_matrix((np.ones(n), (np.arange(n), ids)), shape=(n, n_ids))
J.todense()
matrix([[0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 1., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 1., 0., 0., 0., 0., 0., 0., 0., 0.],
        [1., 0., 0., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 1., 0., 0., 0.],
        [0., 0., 0., 0., 0., 1., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.],
        [0., 0., 1., 0., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 1., 0., 0., 0., 0., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 1., 0., 0.],
        [0., 0., 0., 0., 0., 0., 0., 0., 0., 1.],
        [0., 0., 0., 0., 0., 0., 0., 0., 1., 0.]])
(J.T @ J).todense()
matrix([[10.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  8.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  7.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  0., 13.,  0.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  5.,  0.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  0., 13.,  0.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  0.,  9.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  0.,  0., 15.,  0.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., 14.,  0.],
        [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  6.]])

15.2 Matrix multiplication with diagonal matrices

Doing element-wise multiplication with broadcasting is considerably faster than constructing a diagonal matrix and doing matrix multiplication

n = int(3 * 1e3)
a = rng.normal(size=(n, n))
b = rng.normal(size=n)

Note that to construct a diagonal matrix out of a vector, use the function np.diag. np.diag can also extract the diagonal elements from a matrix and convert them into a vector.

diag_b = np.diag(b)
diag_b
array([[ 0.18434288,  0.        ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.3060658 ,  0.        , ...,  0.        ,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.75518856, ...,  0.        ,
         0.        ,  0.        ],
       ...,
       [ 0.        ,  0.        ,  0.        , ..., -0.07587435,
         0.        ,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
        -0.54171436,  0.        ],
       [ 0.        ,  0.        ,  0.        , ...,  0.        ,
         0.        , -1.09279164]])
%%timeit
diag_b @ a
115 ms ± 5.24 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
%%timeit
b[:, None] * a
8.14 ms ± 328 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%%timeit
a @ diag_b
109 ms ± 819 μs per loop (mean ± std. dev. of 7 runs, 10 loops each)
%%timeit
b[None, :] * a
8.15 ms ± 117 μs per loop (mean ± std. dev. of 7 runs, 100 loops each)