Linear Operators - Scientific Computing with Python

%pylab inline

import time

import numpy as np
import scipy as sp
import scipy.sparse as sparse
import scipy.sparse.linalg as sla
from scipy.sparse.linalg import LinearOperator, aslinearoperator

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

(1)

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

v.sum()

3.055759776070151

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.00027561187744140625 sec.
time for dense: 0.0016298294067382812 sec.
LinearOperator speedup = 5.913494809688581 x
2.1907600843918615e-12

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

8.881784197001252e-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([[2.1794517 , 0.80902747],
       [2.1794517 , 0.80902747],
       [2.1794517 , 0.80902747],
       [2.1794517 , 0.80902747],
       [2.1794517 , 0.80902747]])

np.sum(x, axis=0)

array([2.1794517 , 0.80902747])

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 (computes A @ x)
rmatvec (computes A.T @ x)
matmat (computes A @ B)
rmatmat (computes A.T @ B) as well as the shape of the operator. Note that you don’t need to define all the functions (matvec is 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.float64
)

# 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-25-e8ef22334d5e> in <module>
      1 x = np.random.rand(m)
----> 2 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.05001422, 5.05001422, 5.05001422, 5.05001422, 5.05001422,
       5.05001422, 5.05001422, 5.05001422, 5.05001422, 5.05001422])

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

B2 @ np.random.rand(m)

array([60.71496565, 60.71496565, 60.71496565, 60.71496565, 60.71496565,
       60.71496565, 60.71496565, 60.71496565, 60.71496565, 60.71496565])

type(B1), type(B2)

(scipy.sparse.linalg.interface._SumLinearOperator,
 scipy.sparse.linalg.interface._ProductLinearOperator)

Exercises¶

Define a LinearOperator that “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.float64
)

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

5.087681048627601e-16

Define a LinearOperator that returns the differences in a vector. i.e. A is a n-1 x n operator, where (A @ x)[i] is x[i+1] - x[i]. You may want to look at np.diff.

## Your code here

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