Note
Go to the end to download the full example code
Matrix Multiplication#
This example shows how to use the pylops.MatrixMult operator
to perform Matrix inversion of the following linear system.
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
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.
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()
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()
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()
And we can also use a sparse matrix from the scipy.sparse
family of sparse matrices.
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}")
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.251 seconds)