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 numpy as np
from scipy.sparse import rand
from scipy.sparse.linalg import lsqr

import matplotlib.pyplot as plt
import matplotlib.gridspec as pltgs

import pylops

plt.close('all')
# sphinx_gallery_thumbnail_number = 2

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

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

#a = 1
Aop = pylops.MatrixMult(A, dtype='float64')

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([])

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([])

• 
• 

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()

Out:

/home/docs/checkouts/readthedocs.org/user_builds/pylops/envs/v1.5.0/lib/python3.6/site-packages/numpy/core/numeric.py:538: ComplexWarning: Casting complex values to real discards the imaginary part
  return array(a, dtype, copy=False, order=order)

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()

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('A= %s' % Aop.A.todense())
print('A^-1=', Aop.inv().todense())
print('eigs=', Aop.eigs())
print('x= %s' % x)
print('y= %s' % y)
print('lsqr solution xest= %s' % xest)

Out:

/home/docs/checkouts/readthedocs.org/user_builds/pylops/envs/v1.5.0/lib/python3.6/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  SparseEfficiencyWarning)
A= [[0.23251726 0.93205803 0.59023697 0.89898106 0.33666247]
 [0.         0.32477999 0.97069341 0.46595077 0.80157897]
 [0.49118641 0.         0.89410323 0.34767604 0.        ]
 [0.11482078 0.57378998 0.015679   0.86068497 0.        ]
 [0.48522999 0.         0.         0.88895506 0.        ]]
/home/docs/checkouts/readthedocs.org/user_builds/pylops/envs/v1.5.0/lib/python3.6/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:133: SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix format
  SparseEfficiencyWarning)
/home/docs/checkouts/readthedocs.org/user_builds/pylops/envs/v1.5.0/lib/python3.6/site-packages/scipy/sparse/linalg/dsolve/linsolve.py:203: SparseEfficiencyWarning: spsolve is more efficient when sparse b is in the CSC matrix format
  'is in the CSC matrix format', SparseEfficiencyWarning)
A^-1= [[ 3.48845031 -1.46514609 -0.6274005  -4.83728901  2.16900647]
 [ 2.19028229 -0.91991665 -0.42448569 -1.29437545 -0.31357424]
 [-1.17598641  0.49391327  1.32994116  1.63069147 -1.16861933]
 [-1.90414655  0.79973989  0.34246223  2.64040085 -0.05902097]
 [ 1.64350738  0.55726605 -1.63760684 -2.985122    1.57653162]]
eigs= [ 1.95082375+0.j         -0.22481489+0.71173792j -0.22481489-0.71173792j]
x= [1. 1. 1. 1. 1.]
y= [2.99045578 2.56300314 1.73296567 1.56497474 1.37418506]
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 dims input parameter.

A = np.array([[1., 2.], [4., 5.]])
x = np.array([[1., 1.],
              [2., 2.],
              [3., 3.]])

Aop = pylops.MatrixMult(A, dims=(3,), dtype='float64')
y = Aop*x.flatten()

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

print('A= %s' % A)
print('x= %s' % x)
print('y= %s' % y)
print('lsqr solution xest= %s' % 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.]]