Matrix Multiplication

This example shows how to use the pylops.MatrixMult operator to perform Matrix inversion of the following linear system.

\[\mathbf{y}= \mathbf{A} \mathbf{x}\]

You will see that since this operator is a simple overloading to a numpy.ndarray object, the solution of the linear system can be obtained via both direct inversion (i.e., by means explicit solver such as scipy.linalg.solve or scipy.linalg.lstsq) and iterative solver (i.e., from scipy.sparse.linalg.lsqr).

Note that in case of rectangular \(\mathbf{A}\), an exact inverse does not exist and a least-square solution is computed instead.

import warnings

import matplotlib.gridspec as pltgs
import matplotlib.pyplot as plt
import numpy as np
from scipy.sparse import rand
from scipy.sparse.linalg import lsqr

import pylops

plt.close("all")
warnings.filterwarnings("ignore")
# sphinx_gallery_thumbnail_number = 2

Let’s define the sizes of the matrix \(\mathbf{A}\) (N and M) and fill the matrix with random numbers

N, M = 20, 20
A = np.random.normal(0, 1, (N, M))
Aop = pylops.MatrixMult(A, dtype="float64")

x = np.ones(M)

We can now apply the forward operator to create the data vector \(\mathbf{y}\) and use / to solve the system by means of an explicit solver.

y = Aop * x
xest = Aop / y

Let’s visually plot the system of equations we just solved.

gs = pltgs.GridSpec(1, 6)
fig = plt.figure(figsize=(7, 3))
ax = plt.subplot(gs[0, 0])
ax.imshow(y[:, np.newaxis], cmap="rainbow")
ax.set_title("y", size=20, fontweight="bold")
ax.set_xticks([])
ax.set_yticks(np.arange(N - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 1])
ax.text(
    0.35,
    0.5,
    "=",
    horizontalalignment="center",
    verticalalignment="center",
    size=40,
    fontweight="bold",
)
ax.axis("off")
ax = plt.subplot(gs[0, 2:5])
ax.imshow(Aop.A, cmap="rainbow")
ax.set_title("A", size=20, fontweight="bold")
ax.set_xticks(np.arange(N - 1) + 0.5)
ax.set_yticks(np.arange(M - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 5])
ax.imshow(x[:, np.newaxis], cmap="rainbow")
ax.set_title("x", size=20, fontweight="bold")
ax.set_xticks([])
ax.set_yticks(np.arange(M - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
plt.tight_layout()

gs = pltgs.GridSpec(1, 6)
fig = plt.figure(figsize=(7, 3))
ax = plt.subplot(gs[0, 0])
ax.imshow(x[:, np.newaxis], cmap="rainbow")
ax.set_title("xest", size=20, fontweight="bold")
ax.set_xticks([])
ax.set_yticks(np.arange(M - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 1])
ax.text(
    0.35,
    0.5,
    "=",
    horizontalalignment="center",
    verticalalignment="center",
    size=40,
    fontweight="bold",
)
ax.axis("off")
ax = plt.subplot(gs[0, 2:5])
ax.imshow(Aop.inv(), cmap="rainbow")
ax.set_title(r"A$^{-1}$", size=20, fontweight="bold")
ax.set_xticks(np.arange(N - 1) + 0.5)
ax.set_yticks(np.arange(M - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
ax = plt.subplot(gs[0, 5])
ax.imshow(y[:, np.newaxis], cmap="rainbow")
ax.set_title("y", size=20, fontweight="bold")
ax.set_xticks([])
ax.set_yticks(np.arange(N - 1) + 0.5)
ax.grid(linewidth=3, color="white")
ax.xaxis.set_ticklabels([])
ax.yaxis.set_ticklabels([])
plt.tight_layout()
  • y, A, x
  • xest, A$^{-1}$, y

Let’s also plot the matrix eigenvalues

plt.figure(figsize=(8, 3))
plt.plot(Aop.eigs(), "k", lw=2)
plt.title("Eigenvalues", size=16, fontweight="bold")
plt.xlabel("#eigenvalue")
plt.xlabel("intensity")
plt.tight_layout()
Eigenvalues

We can also repeat the same exercise for a non-square matrix

N, M = 200, 50
A = np.random.normal(0, 1, (N, M))
x = np.ones(M)

Aop = pylops.MatrixMult(A, dtype="float64")
y = Aop * x
yn = y + np.random.normal(0, 0.3, N)

xest = Aop / y
xnest = Aop / yn

plt.figure(figsize=(8, 3))
plt.plot(x, "k", lw=2, label="True")
plt.plot(xest, "--r", lw=2, label="Noise-free")
plt.plot(xnest, "--g", lw=2, label="Noisy")
plt.title("Matrix inversion", size=16, fontweight="bold")
plt.legend()
plt.tight_layout()
Matrix inversion

And we can also use a sparse matrix from the scipy.sparse family of sparse matrices.

N, M = 5, 5
A = rand(N, M, density=0.75)
x = np.ones(M)

Aop = pylops.MatrixMult(A, dtype="float64")
y = Aop * x
xest = Aop / y

print(f"A= {Aop.A.todense()}")
print(f"A^-1= {Aop.inv().todense()}")
print(f"eigs= {Aop.eigs()}")
print(f"x= {x}")
print(f"y= {y}")
print(f"lsqr solution xest= {xest}")

Out:

A= [[0.09004394 0.07465359 0.         0.87688481 0.47238513]
 [0.79938802 0.85640696 0.71661781 0.85258904 0.        ]
 [0.         0.14734131 0.56655652 0.16975416 0.        ]
 [0.29282169 0.46547813 0.90723735 0.62796526 0.        ]
 [0.50711649 0.9563702  0.77100806 0.01124646 0.        ]]
A^-1= [[  0.          13.07319194  29.81875979 -25.7426719   -3.77140228]
 [  0.         -10.01880875 -25.46927387  20.41578795   4.00446478]
 [  0.           3.89143177  12.16350133  -8.55004722  -1.19720943]
 [  0.          -4.29168649 -12.5984235   10.81561744   0.51995029]
 [  2.11691676   7.05799755  21.72748426 -18.39641109  -0.87913918]]
eigs= [ 1.91944139+0.j          0.50019806+0.j         -0.17658613-0.24554532j]
x= [1. 1. 1. 1. 1.]
y= [1.51396747 3.22500183 0.88365199 2.29350244 2.24574121]
lsqr solution xest= [1. 1. 1. 1. 1.]

Finally, in several circumstances the input model \(\mathbf{x}\) may be more naturally arranged as a matrix or a multi-dimensional array and it may be desirable to apply the same matrix to every columns of the model. This can be mathematically expressed as:

\[\begin{split}\mathbf{y} = \begin{bmatrix} \mathbf{A} \quad \mathbf{0} \quad \mathbf{0} \\ \mathbf{0} \quad \mathbf{A} \quad \mathbf{0} \\ \mathbf{0} \quad \mathbf{0} \quad \mathbf{A} \end{bmatrix} \begin{bmatrix} \mathbf{x_1} \\ \mathbf{x_2} \\ \mathbf{x_3} \end{bmatrix}\end{split}\]

and apply it using the same implementation of the pylops.MatrixMult operator by simply telling the operator how we want the model to be organized through the otherdims input parameter.

A = np.array([[1.0, 2.0], [4.0, 5.0]])
x = np.array([[1.0, 1.0], [2.0, 2.0], [3.0, 3.0]])

Aop = pylops.MatrixMult(A, otherdims=(3,), dtype="float64")
y = Aop * x.ravel()

xest, istop, itn, r1norm, r2norm = lsqr(Aop, y, damp=1e-10, iter_lim=10, show=0)[0:5]
xest = xest.reshape(3, 2)

print(f"A= {A}")
print(f"x= {x}")
print(f"y={y}")
print(f"lsqr solution xest= {xest}")

Out:

A= [[1. 2.]
 [4. 5.]]
x= [[1. 1.]
 [2. 2.]
 [3. 3.]]
y=[ 5.  7.  8. 14. 19. 23.]
lsqr solution xest= [[1. 1.]
 [2. 2.]
 [3. 3.]]

Total running time of the script: ( 0 minutes 1.166 seconds)

Gallery generated by Sphinx-Gallery