Linear Operators
Contents
Linear Operators¶
%pylab inline
import scipy as sp
import scipy.sparse as sparse
import scipy.sparse.linalg as sla
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([7.21728315, 7.21728315, 7.21728315, 7.21728315, 7.21728315,
7.21728315, 7.21728315, 7.21728315, 7.21728315, 7.21728315])
v.sum()
7.217283145719205
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.00042247772216796875 sec.
time for dense: 0.016211271286010742 sec.
LinearOperator speedup = 38.371896162528216 x
4.493866839778177e-13
Linear Operators from Functions¶
We can also specify linear operators through functions.
Another way to characterize the action of a matrix containing all 1s is that it sums up all the entries in a vector of length n, and repeats the sum in every entry in a vector of length m. For simplicity, we’ll assume that m = n.
n = 5
a = 1.0
b = 2.0
x = np.random.randn(n)
y = np.random.randn(n)
np.sum(a*x + b*y) - (a*np.sum(x) + b*np.sum(y))
2.220446049250313e-16
# a function that sums the entries in a vector and repeats the sum in every entry of a vector of the same shape.
# works on columns of matrices as well
A = lambda X : np.sum(X, axis=0).reshape(1,-1).repeat(X.shape[0], axis=0)
x = np.random.rand(5, 2)
A(x)
array([[3.23097895, 2.69792801],
[3.23097895, 2.69792801],
[3.23097895, 2.69792801],
[3.23097895, 2.69792801],
[3.23097895, 2.69792801]])
np.sum(x, axis=0)
array([3.23097895, 2.69792801])
The problem is that this function doesn’t play nicely with other linear operators. In order to do that, we wrap the function in the LinearOperator class. For a LinearOperator A, We define the functions
matvec(computesA @ x)rmatvec(computesA.T @ x)matmat(computesA @ B)rmatmat(computesA.T @ B) as well as the shape of the operator. Note that you don’t need to define all the functions (matvecis most important), but you may get errors in certain situations (e.g. taking transpose) if you don’t.
# works on square matrices
Afun = lambda X : np.sum(X, axis=0).reshape(1,-1).repeat(X.shape[0], axis=0)
m = 10 # linear operator of size 10
A = LinearOperator(
shape = (m,m),
matvec = Afun,
rmatvec = Afun,
matmat = Afun,
rmatmat = Afun,
dtype=np.float
)
# the function itself isn't compatible with numpy matrix-vector multiplication
x = np.random.rand(m)
Afun @ x
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
<ipython-input-14-cf1cec2c3dc5> in <module>
1 # the function itself isn't compatible with numpy matrix-vector multiplication
2 x = np.random.rand(m)
----> 3 Afun @ x
ValueError: matmul: Input operand 0 does not have enough dimensions (has 0, gufunc core with signature (n?,k),(k,m?)->(n?,m?) requires 1)
# the linear operator is
A @ x
array([5.06836928, 5.06836928, 5.06836928, 5.06836928, 5.06836928,
5.06836928, 5.06836928, 5.06836928, 5.06836928, 5.06836928])
Composition¶
Because of linearity, sums and products of Linear operators can also have nice properties. This is encoded in the LinearOperator class:
B1 = A + A # acts as the matrix with 2 in every entry
B2 = A @ A # acts that the matrix with m in every entry
B1 @ np.random.rand(m)
array([9.49175069, 9.49175069, 9.49175069, 9.49175069, 9.49175069,
9.49175069, 9.49175069, 9.49175069, 9.49175069, 9.49175069])
B2 @ np.random.rand(m)
array([60.27238783, 60.27238783, 60.27238783, 60.27238783, 60.27238783,
60.27238783, 60.27238783, 60.27238783, 60.27238783, 60.27238783])
type(B1), type(B2)
(scipy.sparse.linalg.interface._SumLinearOperator,
scipy.sparse.linalg.interface._ProductLinearOperator)
Exercises¶
Define a
LinearOperatorthat “centers” a vector:A: x -> x - mean(x). i.e. we subtract the mean of the vector from every entry of the vector.
## Your code here
Afun = lambda x : x - np.mean(x)
m,n = 10,10
A = LinearOperator(
shape = (m,n),
matvec=Afun,
dtype=np.float
)
x = np.random.randn(n)
np.linalg.norm(Afun(x) - (x - np.mean(x)))
np.linalg.norm(A @ x - (x - np.mean(x)))
0.0
# check numerical linearity
x, y = np.random.randn(n), np.random.randn(n)
a, b = 1.5, 2.0
np.linalg.norm(A @ (a*x + b*y) - (a*A@x + b*A@y))
1.1281198993241301e-15
Define a
LinearOperatorthat returns the differences in a vector. i.e.Ais an-1 x noperator, where(A @ x)[i]isx[i+1] - x[i]. You may want to look atnp.diff.
## Your code here
# check numerical linearity
x, y = np.random.randn(n), np.random.randn(n)
a, b = 1.5, 2.0
np.linalg.norm(np.diff(a *x + b * y) - (a*np.diff(x) + b*np.diff(y)))
4.577566798522237e-16
Afun = lambda x : np.diff(x)
A = LinearOperator(
shape = (n-1,n),
matvec = np.diff,
dtype=np.float
)
x = np.random.randn(n)
np.linalg.norm(A @ x - np.diff(x))
0.0
In both the above exercises, you could define these linear operators using either functions or a combination of dense/sparse matrices. If you completed an exercise above using a function, write a second version which uses dense/sparse matrices. If you completed an exercise above using dense/sparse matrices, write a second version that uses functions.
## Your code here