Notebook

In [ ]:

from IPython import display

Sascha Spors, Professorship Signal Theory and Digital Signal Processing, Institute of Communications Engineering (INT), Faculty of Computer Science and Electrical Engineering (IEF), University of Rostock, Germany

Data Driven Audio Signal Processing - A Tutorial with Computational Examples¶

Winter Semester 2023/24 (Master Course #24512)

Feel free to contact lecturer frank.schultz@uni-rostock.de

Exercise 4: Singular Value Decomposition (SVD)¶

Objectives¶

SVD
Four subspaces of a matrix in SVD domain
Rank-1 matrix superposition

Special Python Packages¶

Some convenient functions are found in scipy.linalg, some in numpy.linalg

Highly Recommended Material¶

Kenji Hiranabe's Graphic notes on Gilbert Strang's "Linear Algebra for Everyone"
Gilbert Strang's
Steven L. Brunton's video lecture(s) on SVD

Some Initial Python Stuff¶

In [ ]:

import numpy as np
from scipy.linalg import svd, diagsvd, inv, pinv, null_space, norm
from numpy.linalg import matrix_rank

np.set_printoptions(precision=3, floatmode="fixed", suppress=True)

rng = np.random.default_rng(1234)

# we might convince ourselves that all works for complex data as well
# then the ^H operator (conj().T) needs to be used instead of just .T
use_complex = False

Generate a Matrix with desired rank and dimensions¶

In [ ]:

M = 7  # number of rows
N = 4  # number of cols
# set desired rank, here N-1 for a rank deficient matrix,
# i.e. neither full column nor full row rank
rank = min(M, N) - 1

if False:  # another small toy example
    M = 4  # number of rows
    N = 3  # number of cols
    rank = 2

print("desired rank of A:", rank)

if use_complex:
    dtype = "complex128"
    A = np.zeros([M, N], dtype=dtype)
    for i in range(rank):
        col = rng.integers(-3, 3, M) + 1j * rng.integers(-3, 3, M)
        row = rng.integers(-3, 3, N) + 1j * rng.integers(-3, 3, N)
        A += np.outer(col, row)  # superposition of rank-1 matrices
else:
    dtype = "int"
    A = np.zeros([M, N], dtype=dtype)
    for i in range(rank):
        col = rng.integers(-3, 3, M)
        row = rng.integers(-3, 3, N)
        A += np.outer(col, row)  # superposition of rank-1 matrices
# check if rng produced desired matrix rank
print("        rank of A:", matrix_rank(A))
print("A =\n", A)

We could otherwise consider the matrix A on slide 4/30 from the brilliant A 2020 Vision of Linear Algebra, Slides, 2020. This is the case we have discussed in detail in the 2nd tutorial.

In [ ]:

if True:
    A = np.array([[1, 4, 5], [3, 2, 5], [2, 1, 3]])
    M, N = A.shape[0], A.shape[1]
    rank = matrix_rank(A)
    print("rank r =", rank)
    print("A=\n", A)
    [U, s, Vh] = svd(A)
    V = Vh.T.conj()

    print(
        "column space is spanned by plane in 3D space with these two vectors (1st two columns):"
    )
    print(A[:, :rank])
    print(
        "only two independent columns (the rank tells us this already), because"
    )
    print(
        "A[:,0,None] + A[:,1,None] == A[:,2,None]:\n",
        A[:, 0, None] + A[:, 1, None] == A[:, 2, None],
    )

    print("left null space is spanned by line in 3D space with this vector:")
    lns = np.array([[-1], [7], [-10]])
    print(lns)
    print(
        "check that we get the zero vector for left null space stuff:\n",
        A.T @ lns,
    )

    print("null space is spanned by line in 3D space with this vector:")
    ns = np.array([[-1], [-1], [1]])
    print(ns)
    print("check that we get the zero vector for null space stuff:\n", A @ ns)

    print("row space is spanned by plane in 3D space with two vectors")
    rs = np.array([[1, 0], [0, 1], [1, 1]])
    print("rs=\n", rs)
    print("project rs into V space to get weights for the V vectors")
    print(Vh @ rs)  # third weights are zero, as these belong to nullspace
    print("linear combinations of 1st and 2nd V vectors yield 2 rs vectors:")
    print(V @ (Vh @ rs)[:, 0, None])
    print(V @ (Vh @ rs)[:, 1, None])

    print("we could check that the null space, the left null space,")
    print("the column space and the row space that we came up here")
    print("is also (and actually in a nicer way, because")
    print(
        "SVD orthonormality holds) spanned via the V and U matrices of the SVD"
    )
    print("see code below")

SVD of Matrix A¶

In [ ]:

[U, s, Vh] = svd(A)
S = diagsvd(s, M, N)  # for full SVD the matrix S has same dim as A
V = Vh.conj().T

print("U =\n", U)
# number of non-zero sing values along diag must match rank
print("\nnon-zero singular values: ", s[:rank], "\n")
print("S =\n", S)
print("V =\n", V)

In [ ]:

from IPython.display import Image

# image from the brilliant project Kenji Hiranabe: Graphic notes on
# Gilbert Strang's "Linear Algebra for Everyone"
# found at https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra
# CC0-1.0 License
Image(
    "https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra/raw/2357e987993ba88eb34bbe16e038ce1b150c4878/SVD.png"
)

Four Subspaces in the SVD Domain¶

We denote

column space $C(\mathbf{A})$ , aka image, range
left null space $N(\mathbf{A}^\mathrm{H})$
row space $C(\mathbf{A}^\mathrm{H})$
null space $N(\mathbf{A})$ , aka kernel

The vectors of column space $C(\mathbf{A})$ and left null space $N(\mathbf{A}^\mathrm{H})$ are related to the matrix $\mathbf{U}$ (matrix output).

The vectors of row space $C(\mathbf{A}^\mathrm{H})$ and null space $N(\mathbf{A})$ are related to the matrix $\mathbf{V}$ (matrix input).

Since $\mathbf{U}$ and $\mathbf{V}$ are unitary matrices, they fulfill orthonormality and therefore the properties of orthogonal subspaces

$C(\mathbf{A}) \perp N(\mathbf{A}^\mathrm{H})$

$C(\mathbf{A}^\mathrm{H}) \perp N(\mathbf{A})$

immediately can be deduced. In other words, since we know by the property of SVD, that the matrices $\mathbf{U}$ and $\mathbf{V}$ are orthonormal, we see that

column space $C(\mathbf{A})$ is orthogonal to left null space $N(\mathbf{A}^\mathrm{H})$ , since both spaces are spanned by the dedicated vectors in $\mathbf{U}$
row space $C(\mathbf{A}^\mathrm{H})$ is orthogonal to null space $N(\mathbf{A})$ , since both spaces are spanned by the dedicated vectors in $\mathbf{V}$

In [ ]:

# image from the brilliant project Kenji Hiranabe: Graphic notes on
# Gilbert Strang's "Linear Algebra for Everyone"
# found at https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra
# CC0-1.0 License
Image(
    "https://github.com/kenjihiranabe/The-Art-of-Linear-Algebra/raw/2357e987993ba88eb34bbe16e038ce1b150c4878/4-Subspaces.png"
)

also see the nice graphics of the 4 subspaces with indicated SVD vectors in this wiki http://mlwiki.org/index.php/Four_Fundamental_Subspaces

In [ ]:

# all stuff that is in matrix U (matrix output)
print("U =\n", U)

# column space C(A)
print("\ncolumn space (orthogonal to left null space):")
print(U[:, :rank])

# left null space, if empty only the 0 vector is the left nullspace
print("left null space (orthogonal to column space):")
print(U[:, rank:])

print("###")

# all stuff that is in matrix V (matrix input)
print("\nV =\n", V)

# row space
print("\nrow space (orthogonal to null space):")
print(V[:, :rank])

# null space N(A), if empty only the 0 vector is the null space
print("null space (orthogonal to row space):")
print(V[:, rank:])

In [ ]:

# we use null_space() function

# this very often yields the same vectors as doing the SVD manually:
print("null space: \n", null_space(A))

# this might be different from U[:,rank:], but spans the same space:
print("left null space: \n", null_space(A.conj().T))
# we might want to check this (how should we do this actually?!)

Pure Row Space --> Column Space¶

In [ ]:

# x as linear combination
# of the first two right singular vectors
# i.e. from the row space:
x = 1 * V[:, 0, None] + 2 * V[:, 1, None]
# let matrix A act on x
print(A @ x)
# result must be identical with the linear combination
# of the first two left singular vectors
# i.e. from the column space, weighted by their dedicated singular values and
# same weights 1,2 as above
print(1 * S[0, 0] * U[:, 0, None] + 2 * S[1, 1] * U[:, 1, None])

Row Space + Null Space --> Column Space + Null Vector¶

In [ ]:

# x vector as a linear combination of row space and null space
x_r = 3 * V[:, 0, None]  # from row space
x_n = 4 * V[:, rank, None]  # from null space
print("x from pure null space must yield zero vector in the left nullspace:")
print("A @ x_n =\n", A @ x_n)
x = x_r + x_n  # linear combination
print("let matrix A act on x:")
print(A @ x)
print(
    "note: we transpose the following identical results to have convenient row display"
)
print((A @ x).T)
print((A @ x_r).T)
print((3 * S[0, 0] * U[:, 0, None]).T)

Superposition of Rank-1 Matrices¶

Eckhart-Young theorem, cf. https://en.wikipedia.org/wiki/Low-rank_approximation

In [ ]:

Ar = np.zeros([M, N], dtype="complex128")
for r in range(rank):
    print("\nSum of", r + 1, "rank-1 matrices -> check (A-Ar)")
    Ar += S[r, r] * np.outer(U[:, r], V[:, r].conj())
    print(
        "Frobenius norm (i.e. root(sum squared sing val)): ",
        norm(A - Ar, "fro"),
    )
    print("Spectral norm / L2 norm (i.e. sigma1)", norm(A - Ar, 2))
    print("allclose(A, Ar):", np.allclose(A, Ar))

In the last case Ar is fully reconstructed from rank-1 matrices, yielding Frobenius and spectral norm for A-Ar of 0 (numerical precision does not give us zero). We can also check with allclose()

Copyright¶

the notebooks are provided as Open Educational Resources
the text is licensed under Creative Commons Attribution 4.0
the code of the IPython examples is licensed under the MIT license
feel free to use the notebooks for your own purposes
please attribute the work as follows: Frank Schultz, Data Driven Audio Signal Processing - A Tutorial Featuring Computational Examples, University of Rostock ideally with relevant file(s), github URL https://github.com/spatialaudio/data-driven-audio-signal-processing-exercise, commit number and/or version tag, year.