The R package numpy.linalg contains the function svd to determine numerically the singular value decomposition.
First we recalculate the example of the learning material:

>>> import numpy as np
>>> import numpy.linalg as la
>>> A = np.array([[2,0],[0,-3]]); print(A)
[[ 2  0]
 [ 0 -3]]
>>> U, S, Vt = la.svd(A, full_matrices=True) # singular value decomposition (SVD) of A
>>> print(U)           # matrix U in SVD
[[0. 1.]
 [1. 0.]]
>>> print(np.diag(S))  # diagonal matrix with singular values 
[[3. 0.]
 [0. 2.]]
>>> print(Vt)          # transposed of matrix V in SVD
[[-0. -1.]
 [ 1. 0.]]
>>> print(U @ np.diag(S) @ Vt - A)  # check of svd(A)
[[0. 0.]
 [0. 0.]]

The following example of a singular matrix illustrates the difference between an exact and numerical computation. \[\begin{aligned}\matrix{1 & -1\\ -1 & 1} &= \matrix{-\frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2} \\ \frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2}}\matrix{2 & 0\\ 0 & 0}\matrix{-\frac{1}{2}\sqrt{2} & \frac{1}{2}\sqrt{2} \\ -\frac{1}{2}\sqrt{2} & -\frac{1}{2}\sqrt{2}}\\ \\ &=\frac{1}{2}\sqrt{2} \matrix{-1 & 1\\ 1 & 1}\matrix{2 & 0\\ 0 & 0}\matrix{-1 & 1\\ -1 & -1}\frac{1}{2}\sqrt{2} \end{aligned} \]

>>> B = np.array([[1, -1],[-1,1]]); print(B)
[[ 1 -1]
 [-1  1]]
>>> U, S, Vt = la.svd(B, full_matrices=True) # singuler value decomposition
>>> print(U)                       # matrix U in SVD
[[-0.70710678 0.70710678]
 [ 0.70710678 0.70710678]]
>>> print(np.diag(S))              # diagonal matrix with singular values 
[[2.00000000e+00 0.00000000e+00]
 [0.00000000e+00 3.35470445e-17]]
>>> print(Vt)                      # transposed of matrix V in SVD
[[-0.70710678  0.70710678]
 [-0.70710678 -0.70710678]]
>>> print(U @ np.diag(S) @ Vt - B) # check of svd(B)
[[-3.33066907e-16  4.44089210e-16]
 [ 0.00000000e+00 -2.22044605e-16]]
>>> print(np.round(U @ np.diag(S) @ Vt - B, 3)) # rounding
[[-0.  0.]
 [ 0. -0.]]

Next comes an example of digital image compression, using a pattern that consists of black and white squares in a \(9\times 9\) grid. We specify the figure with a matrix consisting of zeros (white squares) and ones (for black boxes: \[M = \matrix{ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0}\] In R we can visualise this pattern. We specify the figure with a matrix consisting of zeros (white squares) and ones (for black boxes: \[M = \matrix{ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 1 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 1\\ 0 & 1 & 0 & 1 & 1 &1 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 1 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0}\] In MATLAB we can visualise this pattern.

>>> import numpy.linalg as la
>>> import matplotlib.pyplot as plt
>>> M = np.array(
  [[0, 0, 0, 0, 1, 0, 0, 0, 0],
  [0, 0, 0, 1, 0, 1, 0, 0, 0],
  [0, 0, 1, 0, 0, 0, 1, 0, 0],
  [0, 1, 0, 1, 1, 1, 0, 1, 0],
  [1, 0, 0, 1, 0, 1, 0, 0, 1],
  [0, 1, 0, 1, 1, 1, 0, 1, 0],
  [0, 0, 1, 0, 0, 0, 1, 0, 0],
  [0, 0, 0, 1, 0, 1, 0, 0, 0],
  [0, 0, 0, 0, 1, 0, 0, 0, 0]])  #  patroon
>>> plt.imshow(np.ones_like(M) - M, cmap='gray')
>>> plt.show()

Suppose we want to send someone this pattern. Then we can send the matrix \(M\), that is. 81 numbers one and zero. But you can send fewer numbers when you use the singular value decomposition, while the image can still be reconstructed.

The eigenvalues of the matrix \(M\) are as follows in increasing order:

>>> print(np.flip(np.round(np.sort(vals), 3))
 3.709 2. 1.194 0. 0. -0. -0. -0.903 -2. ])

Thus, the singular values of \(M\) are in descending order, rounded to 1 decimal place, equal to \(3.7\), \(2.0\), \(1.2\), \(0.9\) and \(0\).

We approximate \(M\) now with the first terms of the singular value decomposition. But instead of using the numerical values of the decomposition for image display we round the elements of the matrix elements to whole numbers \(0\) and \(1\) to reconstruct the pattern.

Suppose we only use the highest singular value \(3.7\). Then the picture looks already a bit like the original figure:

>>> U, S, Vt = la.svd(M, full_matrices=True) # singuler value decomposition (SVD) of M
>>> S1 = np.zeros((9,9)); S1[0,0] = S[0] 
>>> M1 = np.round( U @ S1 @ Vt)
>>> plt.imshow(np.ones_like(M1) - M1, cmap='gray')
>>> plt.show()

To do so we only need to transmit \(1+2\times 9 = 19\) numeric values, namely the singular value and the two vectors \(\vec{u}_1\) and \(\vec{v}_1\).
If we also include the singular value \(2\), we get even back the original figure.

>>> S3 = np.zeros((9,9)); S3[0,0] = S[0]; S3[1,1] = S[1]; S3[2,2] = S[2]
>>> M3 = np.round(U @ S3 @ Vt)
>>> plt.imshow(np.ones_like(M3) - M3, cmap='gray')
>>> plt.show()

Now it suffices to transmit \(3+6\times 9 = 57\) numerical values, namely \(3\) eigenvalues and \(2\times 3= 6\) vectors. This does not seem so much less than \(81\), it is already a reduction of approximately \(30\%\)

\(\phantom{xx}\)

The gain is more spectacular when transmitting large digital images, for example in the picture commonly used in computer vision from \(512\times 512\) pixels.

We read the image into the software environment, calculate the singular value decomposition of the corresponding \(512\times 512\) matrix with grayscale values, and approximate the original image with 8, 16, 32, 64, and 128 terms from the singular value decomposition.

>>> import numpy as np
>>> import numpy.linalg as la
>>> import matplotlib.pyplot as plt
>>> from imageio import imread
>>> img = imread('D:/lena512.bmp')  # load the picture
>>> plt.imshow(img, cmap='gray')
>>> plt.show()    # display the above photo (with axes)
>>> [U,S,Vt] = la.svd(img)  # singulierewaardenontbinding
>>> img8 = U[:,0:8] @ np.diag(S[0:8]) @ Vt[0:8,:]         # image compression k=  8
>>> img16 = U[:,0:16] @ np.diag(S[0:16]) @ Vt[0:16,:]     # image compression k= 16
>>> img32 = U[:,0:32] @ np.diag(S[0:32]) @ Vt[0:32,:]     # image compression k= 32
>>> img64 = U[:,0:64] @ np.diag(S[0:64]) @ Vt[0:64,:]     # image compression k= 64
>>> img128 = U[:,0:128] @ np.diag(S[0:128]) @ Vt[0:128,:] # image compression k=128
>>> fig = plt.figure()  # display images in a 3x2 grid
>>> fig.add_subplot(3,2,1)   
>>> plt.axis('off')
>>> plt.imshow(img8, cmap='gray')
>>> fig.add_subplot(3,2,2)
>>> plt.axis('off')
>>> plt.imshow(img16, cmap='gray')
>>> fig.add_subplot(3,2,3)
>>> plt.axis('off')
>>> plt.imshow(img32, cmap='gray')
>>> fig.add_subplot(3,2,4)
>>> plt.axis('off')
>>> plt.imshow(img64, cmap='gray')
>>> fig.add_subplot(3,2,5)
>>> plt.axis('off')
>>> plt.imshow(img128, cmap='gray')
>>> fig.add_subplot(3,2,6)
>>> plt.axis('off')
>>> plt.imshow(img, cmap='gray')
>>> plt.show()

Below we present the images for different values of the number of terms \(k\) extracted from the SVD. Clearly less information you need while still having a picture that differs on first sight little form the original (at \(k=64\) the compression factor is about \(\frac{1}{4}\)). With colour pictures, the gain is even greater.

\(k=8\)	\(k=16\)
\(k=32\)	\(k=64\)
\(k=128\)	Original: \(k=512\)