idss_lecture_02_vectorsmatrices_20232024a.pdf

numpy In : # standard imports import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt %matplotlib inline from jhwutils.matrices import show_matrix_effect, print_matrix plt.rc('figure', figsize=(6.0, 5.0), dpi=240) ...

numpy In : # standard imports import numpy as np import matplotlib as mpl import matplotlib.pyplot as plt %matplotlib inline from jhwutils.matrices import show_matrix_effect, print_matrix plt.rc('figure', figsize=(6.0, 5.0), dpi=240)  salamander salamander axolotl waterdog king woman man 5 * [1,0,0] + 7 * [0,1,0] + 3 * [0,0,1] In : # two 3D vectors (3 element ordered tuples) x = np.array([0,1,1]) y = np.array([4,5,6]) a = 2 print_matrix("a", a) print_matrix("x", x) print_matrix("y", y) x [[0 1 1]] y [[4 5 6]] In : print_matrix("ax", a*x) print_matrix("ay", a*y) print_matrix("x+y", x+y) print_matrix("\|x\|_2", np.linalg.norm(x)) # norm print_matrix("x\cdot y", np.dot(x,y)) # inner product ax [[0 2 2]] ay [[ 8 10 12]] x+y [[4 6 7]] heart_rate systolic diastolic vo2 67 110 72 98 65 111 70 98 64 110 69 97.. [[ -0. 15. 0. ] [ 16. -0.5 32.5] [-16. -0.5 32.5] [ 16. -15. -32.5] [-16. -15. -32.5] [-44. 10. -32.5] [-60. -3. -13. ] [-65. -3. -32.5] [ 44. 10. -32.5] [ 60. -3. -13. ] [ 65. -3. -32.5] [ -0. 15. -32.5]] # some GLSL vec3 pos = vec3(1.0,2.0,0.0); vec3 vel = vec3(0.1,0.0,0.0); pos = pos + vel; In : v1 = np.array([2.0, 1.0]) v2 = np.array([-1.0, 0.5]) v3 = np.array([-0.5, 0.0]) origin = np.array([0,0]) def show_vector(ax, start, end, label='', color=None, **kwargs): vec = np.stack([start, end]) lines = ax.plot(end, end, 'o', label=label, color=color, **kwargs) if color is None: color = lines.get_color() ax.arrow(start, start, end-start, end-start, head_width=0.1, width=0.02, overhang=0.2, length_includes_head=True, color=color, **kwargs) fig = plt.figure() ax = fig.add_subplot(1,1,1) # show the original vectors show_vector(ax, origin, v1, 'V1') show_vector(ax, origin, v2, 'V2') show_vector(ax, origin, v3, 'V3') # show some sums of vectors show_vector(ax, v1, v1+v2, 'V1+V2') show_vector(ax, v1+v2, v1+v2+v3, 'V1+V2+V3') show_vector(ax, v1+v2+v3, v1+v2+v3-v2, 'V1+V2+V3-V2') ax.set_frame_on(False) ax.set_xlim(-3,3) ax.set_ylim(-3,3) ax.set_aspect(1.0) ax.legend() Out: In : fig = plt.figure() ax = fig.add_subplot(1, 1, 1) # show the original vectors show_vector(ax, origin, v1, 'V1') show_vector(ax, origin, v2, 'V2') # show a range of weighted vectors in between alphas = np.linspace(0, 1, 15) for alpha in alphas:  show_vector(ax, origin, alpha * v1 + (1 - alpha) * v2, color='k', alpha=0.25, label='') ax.set_xlim(-2, 3) ax.set_ylim(-2, 3) ax.set_title("Linear interpolation between two 2D vectors") ax.legend() ax.set_aspect(1.0)  np.linalg.norm() In : x = np.array([1.0, 10.0, -5.0]) y = np.array([1.0, -4.0, 8.0]) print_matrix("x", x) print_matrix("y", y) print_matrix("\|x\|", np.linalg.norm(x)) print_matrix("\|y\|", np.linalg.norm(y)) x [[ 1. 10. -5.]] y [[ 1. -4. 8.]] norm In : test_vector = np.array([1, 0, 2, 0.1, -4, 0]) print_matrix("x", test_vector) for norm in [1, 2, np.inf, -np.inf, 0, 0.5, 3, -1]: print_matrix("\|x\|_{{{norm}}}".format(norm=norm), np.linalg.norm(test_vector,norm)) x [[ 1. 0. 2. 0.1 -4. 0. ]] In : x = np.random.normal(0,5,(5,)) # a random vector x_norm = x / np.linalg.norm(x) # a random unit vector  print_matrix("x", x) print_matrix("x_n", x_norm) x [[ 4.64 1.67 -1.37 -2.45 -3.74]] x_n [[ 0.68 0.25 -0.2 -0.36 -0.55]] In : def connect_plot(ax, a, b): for a1,b1 in zip(a,b): #ax.plot([a1, b1],[a1, b1], 'k-', lw=0.25) ax.arrow(a1, a1, b1-a1, b1-a1, head_width=0.05, length_includes_head=True, facecolor='k', zorder=10) In : # show that 2D unit vectors lie on the unit circle x = np.random.normal(0,1,(100,2)) # 100 2D vectors unit_x = (x.T / np.linalg.norm(x, axis=1)).T # a random unit vector # plot the results fig = plt.figure() ax = fig.add_subplot(1, 1, 1) ax.plot(x[:,0], x[:,1], '.', label="Original") ax.plot(unit_x[:,0], unit_x[:,1], 'o', label="Normalised") ax.legend() connect_plot(ax, x, unit_x) ax.set_aspect(1.0) ax.set_frame_on(False) ax.set_title("Euclidean normalisation of vectors") Text(0.5, 1.0, 'Euclidean normalisation of vectors') Out: In : # in a different norm x = np.random.normal(0,1,(100,2)) # 100 2D vectors unit_x = (x.T / np.linalg.norm(x, axis=1, ord=np.inf)).T # a random L_inf unit vector fig = plt.figure() ax = fig.add_subplot(1, 1, 1) connect_plot(ax, x, unit_x) ax.set_title("$L_\infty$ normalisation of vectors") ax.plot(x[:,0], x[:,1], '.', label="Original") ax.set_frame_on(False) ax.plot(unit_x[:,0], unit_x[:,1], 'o', label="Normalised") ax.legend() ax.set_aspect(1.0) In : x = np.array([1, 2, 3, 4]) y = np.array([4, 0, 1, 4]) print_matrix("x", x) print_matrix("y", y) print_matrix("x\cdot y", np.inner(x, y)) # inner product is same as dot product for vectors x [[1 2 3 4]] y [[4 0 1 4]] In : print(np.inner(x,[1,1,1])) # inner product is not defined for vectors of differing dimension --------------------------------------------------------------------------- ValueError Traceback (most recent call last) Cell In , line 1 ----> 1 print(np.inner(x,[1,1,1])) File :180, in inner(*args, **kwargs) ValueError: shapes (4,) and (3,) not aligned: 4 (dim 0) != 3 (dim 0) np.mean(x, axis=0) In : x = np.random.normal(2,1, (3,4)) # 3 rows, 4 columns print_matrix("x - ", x) mu = np.mean(x, axis=0) # the mean vector print_matrix("np.mean(x, axis=0) - ", mu) # verify the computation is the same as doing it "by hand" print_matrix("np.sum(x, axis=0)/x.shape - ", np.sum(x, axis=0)/x.shape)  x - [[2.07 3.08 1.11 1.54] [2.44 2.98 3.97 2.3 ] [3.37 2.15 2.14 1.85]] np.mean(x, axis=0) - [[2.63 2.74 2.4 1.9 ]] np.sum(x, axis=0)/x.shape - [[2.63 2.74 2.4 1.9 ]] In : x_center = x - mu print_matrix("x_center", x_center) mu_c = np.mean(x_center, axis=0) # verify the mean is now all zeros print_matrix("\mu_c", mu_c) x_center [[-0.56 0.35 -1.3 -0.36] [-0.19 0.24 1.56 0.4 ] [ 0.75 -0.59 -0.27 -0.04]] \mu_c [[ 0. 0. 0. -0.]] In : # show the effect of centering a collection of vectors fig = plt.figure() ax = fig.add_subplot(1,1,1) # original ax.scatter(x[:,0], x[:,1], c='C0', label="Original", s=10) # original mean ax.scatter(mu, mu, c='C2', label="Mean", s=40) # centered ax.scatter(x_center[:,0], x_center[:,1], c='C1', label="Centered", s=10) ax.plot([0, mu], [0, mu], c='C2', label='Shift') # draw origin and fix axes ax.set_xlim(-4,4) ax.set_ylim(-4,4) ax.set_frame_on(False) ax.axhline(0, c='k', ls=':') ax.axvline(0, c='k', ls=':') ax.set_aspect(1.0) ax.legend() Out: In : time = np.linspace(0,120,10000) season = np.sin(2*np.pi*time/12.0) temps = np.random.normal(14.0, 4.0, time.shape) + season * 5 fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.plot(time, temps) ax.set_xlabel("Time (months)") ax.set_ylabel("Temp (deg. C)") ax.set_frame_on(False) ax.set_title("Simulated temperature") Text(0.5, 1.0, 'Simulated temperature') Out: In : fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.hist(temps, bins=20) ax.set_frame_on(False) ax.set_title("Histogram of temperatures") ax.set_xlabel("Temp (deg. C)") ax.set_ylabel("Count") Text(0, 0.5, 'Count') Out: In : humidity = np.random.normal(np.tanh(temps*0.15)*60+30, 5, time.shape) fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.set_xlabel("Time (months)") ax.set_ylabel("Relative humidity (%)") ax.set_title("Simulated humidity") ax.set_frame_on(False) ax.plot(time, humidity) [] Out: In : fig = plt.figure() ax = fig.add_subplot(1,1,1) bar = plt.hist2d(temps, humidity, bins=(20,20)); ax.set_xlabel("Temp (deg. C)") ax.set_ylabel("Relative humidity (%)") fig.colorbar(bar[-1]) ax.set_title("2D histogram of temperature and humidity") Text(0.5, 1.0, '2D histogram of temperature and humidity') Out: In : import scipy.special # for the gamma function def sphere_volume(n): return 0.5**n * np.pi**(n/2.0) / scipy.special.gamma(n/2.0+1) # this one is easy... def cube_volume(n): return 1.0 In : x = np.arange(1,50) fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.plot(x, [sphere_volume(xi) for xi in x], 'o-') ax.set_xlabel("Dimension") ax.set_ylabel("Volume") ax.set_frame_on(False) ax.set_title("Volume of sphere as fraction of circumscibed cube vs. dimension") Text(0.5, 1.0, 'Volume of sphere as fraction of circumscibed cube vs. dimension') Out: In : fig = plt.figure() ax = fig.add_subplot(1,1,1) ax.semilogy(x, [sphere_volume(xi) for xi in x]) ax.set_xlabel("Dimension") ax.set_ylabel("Volume") ax.set_frame_on(False) ax.set_title("Volume of sphere as fraction of circumscibed cube vs. dimension (log scale)") Text(0.5, 1.0, 'Volume of sphere as fraction of circumscibed cube vs. dimension (log scale)') Out: In : # let's do an experiment d = 1 # dimensions n = 100_000 # number of points # cube, at origin, n-dimensional points_in_box = np.random.uniform(-1, 1, (n,d)) # select all points with radius < 1 (inside sphere in that box) points_in_sphere = points_in_box[ np.linalg.norm(points_in_box,axis=1) 4 print_matrix("Cx", C @ x) ValueError: matmul: Input operand 1 has a mismatch in its core dimension 0, with gufunc signa ture (n?,k),(k,m?)->(n?,m?) (size 3 is different from 2) In : def matmul(a, b): p, q_a = a.shape q_b, r = b.shape # we can only multiply two matrices if A is p x q and B in q x r assert q_a == q_b # the result is a matrix of size p x r c = np.zeros((p, r)) for i in range(p): for j in range(r): # Note that this can be seen as a simple *weighted sum* # the sum of the ith row of A weighted by the jth column of B c[i, j] = np.sum(a[i, :] * b[:, j]) return c In : a = np.array([[1, 2, -3]]) b = np.array([[1, -1, 1], [2, -2, 2], [3, -3, 3]]) print_matrix("{\\bf A}", a) print_matrix("{\\bf B}", b) c = matmul(a,b) print_matrix("{\\bf A}B", c) {\bf A} [[ 1 2 -3]] {\bf B} [[ 1 -1 1] [ 2 -2 2] [ 3 -3 3]] {\bf A}B [[-4. 4. -4.]] np.dot(a,b) a @ b a @ b In : # verify that this is the same as the built in matrix multiply c_numpy = np.dot(a, b) print_matrix("C_{\\text numpy}", c_numpy) print(np.allclose(c, c_numpy)) c_at = a @ b print_matrix("C_{\\text a @ b}", a @ b) print(np.allclose(c_at, c)) C_{\text numpy} [[-4 4 -4]] True C_{\text a @ b} [[-4 4 -4]] True In : print_matrix("A", A) print_matrix("x", x) print_matrix("Ax", A @ x) # force to a 2D array of m x 1 size print_matrix("Ax", A @ x.reshape(-1, 1)) # note this will still be a 2D array in NumPy # and will appear differently when printed. However, it # has the same *semantics* as a 1D array A [[ 1 -1 1] [ 1 0 0] [ 0 1 1]] x [[1 2 3]] Ax [[2 1 5]] Ax [ ] A.T In : ### Transpose A = np.array([[2, -5], [1, 0], [3, 3]]) print_matrix("A", A) print_matrix("A^T", A.T) A [[ 2 -5] [ 1 0] [ 3 3]] A^T [[ 2 1 3] [-5 0 3]] In : x = np.array([[1,2,3]]) y = np.array([[4,5,6]]) print_matrix("{\\bf x}", x) print_matrix("{\\bf y}", y) {\bf x} [[1 2 3]] {\bf y} [[4 5 6]] In : print_matrix("{\\bf x} \otimes {\\bf y}", np.outer(x,y)) print_matrix("{\\bf x}^T{\\bf y}", x.T @ y) print_matrix("{\\bf x} \\bullet {\\bf y}", np.inner(x,y)) print_matrix("{\\bf x}{\\bf y}^T", x @ y.T) {\bf x} \otimes {\bf y} [[ 4 5 6] [ 8 10 12] [12 15 18]] {\bf x}^T{\bf y} [[ 4 5 6] [ 8 10 12] [12 15 18]] {\bf x} \bullet {\bf y} [] {\bf x}{\bf y}^T [] In : ## Rotation d30 = np.radians(30)  cs = np.cos(d30) ss = np.sin(d30) rot30 = np.array([[cs, ss], [-ss, cs]]) ## Scaling scale_x = np.array([[1,0], [0, 0.5]]) In : show_matrix_effect(rot30) [[ 0.87 0.5 ] [-0.5 0.87]] In : show_matrix_effect(scale_x) [[1. 0. ] [0. 0.5]] In : A = scale_x @ rot30 # product of scale and rotate print(A) show_matrix_effect(A, "Rotate then scale") # rotate, then scale [[ 0.8660254 0.5 ] [-0.25 0.4330127]] Rotate then scale [[ 0.87 0.5 ] [-0.25 0.43]] In : B = rot30 @ scale_x print(B) show_matrix_effect(B, "Scale then rotate") # scale, then rotate # note: Not the same! [[ 0.8660254 0.25 ] [-0.5 0.4330127]] Scale then rotate [[ 0.87 0.25] [-0.5 0.43]] In : x = np.random.normal(0,1,(500, 5)) mu = np.mean(x, axis=0) sigma_cross = ((x - mu).T @ (x - mu)) / (x.shape-1) np.set_printoptions(suppress=True, precision=2) print_matrix("\Sigma_{\\text{cross}}", sigma_cross) \Sigma_{\text{cross}} [[ 1.05 0.09 0.01 -0. -0.07] [ 0.09 0.94 0.01 -0.02 -0. ] [ 0.01 0.01 0.98 -0.03 0.1 ] [-0. -0.02 -0.03 0.93 -0.02] [-0.07 -0. 0.1 -0.02 1.1 ]] np.cov(x) In : # verify it is close to the function provided by NumPy sigma_np = np.cov(x, rowvar=False) print_matrix("\Sigma_{\\text{numpy}}", sigma_np) \Sigma_{\text{numpy}} [[ 1.05 0.09 0.01 -0. -0.07] [ 0.09 0.94 0.01 -0.02 -0. ] [ 0.01 0.01 0.98 -0.03 0.1 ] [-0. -0.02 -0.03 0.93 -0.02] [-0.07 -0. 0.1 -0.02 1.1 ]] In : # These are helper functions. You do not need to undersatnd these in detail. from matplotlib.patches import Ellipse def eigsorted(cov): vals, vecs = np.linalg.eigh(cov) order = vals.argsort()[::-1] return vals[order], vecs[:, order] def cov_ellipse(ax, x, nstd=1, **kwargs): cov = np.cov(x.T) mu = np.mean(x, axis=0) vals, vecs = eigsorted(cov) theta = np.degrees(np.arctan2(*vecs[:, 0][::-1])) w, h = 2 * nstd * np.sqrt(vals) ell = Ellipse(xy=(mu, mu), width=w, height=h, angle=theta, **kwargs) ax.add_artist(ell) return mu,cov In : fig = plt.figure() ax = fig.add_subplot(1,1,1) x = np.random.normal(0,1,(200,2)) @ np.array([[0.1, 0.5], [-0.9, 1.0]]) ax.scatter(x[:,0], x[:,1], c='C0', label="Original", s=10) cov_ellipse(ax, x[:,0:2], 1, facecolor='none', edgecolor='k') cov_ellipse(ax, x[:,0:2], 2, facecolor='none', edgecolor='k') mu, cov = cov_ellipse(ax, x[:,0:2], 3, facecolor='none', edgecolor='k') ax.set_xlim(-4,4) ax.set_ylim(-4,4) ax.axhline(0) ax.axvline(0) ax.set_frame_on(False) ax.set_aspect(1.0) ax.legend() print_matrix("\mu",mu) print_matrix("\Sigma",cov) \mu [[-0.01 -0.05]] \Sigma [[ 0.79 -0.79] [-0.79 1.19]]  np.diag(x) x In : print_matrix("\\text{diag}(x)", np.diag([1,2,3,3,2,1])) \text{diag}(x) [[1 0 0 0 0 0] [0 2 0 0 0 0] [0 0 3 0 0 0] [0 0 0 3 0 0] [0 0 0 0 2 0] [0 0 0 0 0 1]] In : # an anti-diagonal array is the diagonal of the flipped array print_matrix("\\text{anti-diagonal}", np.fliplr(np.diag([1,2,3]))) \text{anti-diagonal} [[0 0 1] [0 2 0] [3 0 0]]  np.eye(n) In : print_matrix("I", np.eye(3)) I [[1. 0. 0.] [0. 1. 0.] [0. 0. 1.]] In : # your identity never changes anything print_matrix("Ix", np.eye(3) @ np.array([4,5,6]) ) Ix [[4. 5. 6.]] In : print_matrix("AI", np.array([[1,2,3],[4,5,6],[7,8,9]]) @ np.eye(3)) print_matrix("IA", np.eye(3) @ np.array([[1,2,3],[4,5,6],[7,8,9]]) ) AI [[1. 2. 3.] [4. 5. 6.] [7. 8. 9.]] IA [[1. 2. 3.] [4. 5. 6.] [7. 8. 9.]] In : a = 0.5 x = np.array([4,5,6]) aI = np.eye(3) * a print("c=",a) print_matrix("cI", aI) print_matrix("(cI){\\bf x}\n", aI @ x) # the same thing: print_matrix("c{\\bf x}\n",a*x) c= 0.5 cI [[0.5 0. 0. ] [0. 0.5 0. ] [0. 0. 0.5]] (cI){\bf x} [[2. 2.5 3. ]] c{\bf x} [[2. 2.5 3. ]] In : z = np.zeros((4,4)) print(z) [[0. 0. 0. 0.] [0. 0. 0. 0.] [0. 0. 0. 0.] [0. 0. 0. 0.]] In : x = np.array([1,2,3,4]) y = np.array([[1,-1,1], [1,1,-1], [1,1,1], [-1,-1,-1]]) print_matrix("x", x) print_matrix("y", y) print_matrix("0x", z @ x) print() print_matrix("0y", z @ y) x [[1 2 3 4]] y [[ 1 -1 1] [ 1 1 -1] [ 1 1 1] [-1 -1 -1]] 0x [[0. 0. 0. 0.]] 0y [[0. 0. 0.] [0. 0. 0.] [0. 0. 0.] [0. 0. 0.]] In : # tri generates an all ones lower triangular matrix upper = np.tri(4) print_matrix("T_u", upper) # transpose changes a lower triangular to an upper triangular lower = np.tri(4).T print_matrix("T_l", lower) T_u [[1. 0. 0. 0.] [1. 1. 0. 0.] [1. 1. 1. 0.] [1. 1. 1. 1.]] T_l [[1. 1. 1. 1.] [0. 1. 1. 1.] [0. 0. 1. 1.] [0. 0. 0. 1.]]