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.
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.
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.357 seconds)