Linear Operators

Linear Operators#

%pylab inline
import scipy as sp
import scipy.sparse as sparse
import scipy.sparse.linalg as sla

%pylab is deprecated, use %matplotlib inline and import the required libraries.
Populating the interactive namespace from numpy and matplotlib

In linear algebra, a linear transformation, linear operator, or linear map, is a map of vector spaces $T:V \to W$ where $$T(\alpha v_1 + \beta v_2) = \alpha T v_1 + \beta T v_2$$

If you choose bases for the vector spaces $V$ and $W$, you can represent $T$ using a (dense) matrix. However, there are many situations where we may want to represent $T$ in some other format which will allow us to do faster matrix-vector and matrix-matrix multiplications.

The case of a sparse matrix is handled by special matrix formats, but there are also situations in which dense matrices can also be applied quickly.

Low Rank Matrices#

The easiest situation to describe is a low-rank matrix. We saw an example of this when we first saw object-oriented programming in the LinearMap class. SciPy provides a very similar class which can be used to construct arbitrary linear operators, which is called LinearOperator.

The LinearOperator class can be found in scipy.sparse.linalg. The aslinearoperator function lets us to treat dense and sparse arrays as LinearOperators.

from scipy.sparse.linalg import LinearOperator, aslinearoperator

As an example, let’s construct a LinearOperator that acts as the matrix of all ones. This matrix is rank-1 and can be written as $11^T$, where $1$ is a vector of the appropriate dimension.

n = 20
m = 10
onesn = aslinearoperator(np.ones((n,1)))
onesm = aslinearoperator(np.ones((m,1)))
A = onesm @ onesn.T
A

<10x20 _ProductLinearOperator with dtype=float64>

v = np.random.randn(n)
A @ v

array([-2.05634486, -2.05634486, -2.05634486, -2.05634486, -2.05634486,
       -2.05634486, -2.05634486, -2.05634486, -2.05634486, -2.05634486])

v.sum()

-2.0563448646523903

A @ np.eye(A.shape[1])

array([[1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.],
       [1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,
        1., 1., 1., 1.]])

Let’s do a timing comparison between the LinearOperator and a dense matrix.

n = 2000
m = 1000
onesn = aslinearoperator(np.ones((n,1)))
onesm = aslinearoperator(np.ones((m,1)))
Alo = onesm @ onesn.T
Ad = np.ones((m,n)) # dense array of all ones

import time

x = np.random.randn(n)
y = np.empty(m)

t0 = time.time()
y1 = Alo @ x
t1 = time.time()
print("time for linear operator: {} sec.".format(t1 - t0))
tlo = t1 - t0

t0 = time.time()
y2 = Ad @ x
t1 = time.time()
print("time for dense: {} sec.".format(t1 - t0))
td = t1 - t0

print("LinearOperator speedup = {} x".format(td / tlo))
print(np.linalg.norm(y2 - y1))

time for linear operator: 0.00020122528076171875 sec.
time for dense: 0.0075495243072509766 sec.
LinearOperator speedup = 37.51777251184834 x
2.2469334198890886e-13