In [2]:
from __future__ import division
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as pt
import random
In [3]:
def make_random_sparse_matrix(n, row_fill):
nentries = (n*row_fill) // 2 # because of symmetry
data = np.random.randn(nentries)
rows = np.random.random_integers(0, n-1, nentries)
cols = np.random.random_integers(0, n-1, nentries)
import scipy.sparse as sps
coo = sps.coo_matrix((data, (rows, cols)), shape=(n,n))
# NOTE: Cuthill-McKee applies only to symmetric matrices!
return (100*np.eye(n) + np.array(coo.todense() + coo.todense().T))
In [47]:
prec = 1e-10
np.random.seed(15)
random.seed(15)
A = make_random_sparse_matrix(200, 3)
print "%d non-zeros" % len(np.where(np.abs(A)>prec)[0])
pt.figure()
pt.spy(A, marker=",", precision=prec)
Out[47]:
In [48]:
Ainv = la.inv(A)
print "%d non-zeros" % len(np.where(np.abs(Ainv) > prec)[0])
pt.spy(Ainv, marker=",", precision=prec)
Out[48]:
In [49]:
L = la.cholesky(A)
print "%d non-zeros" % len(np.where(np.abs(L) > prec)[0])
pt.spy(L, marker=",", precision=prec)
Out[49]:
Cholesky is often less bad, but in principle affected the same way.
Reducing the fill-in
In [50]:
def degree(mat, row):
return len(np.where(mat[row])[0])
print degree(A, 3)
print degree(A, 4)
print degree(A, 5)
In [51]:
def argmin2(iterable, return_value=False):
it = iter(iterable)
try:
current_argmin, current_min = it.next()
except StopIteration:
raise ValueError("argmin of empty iterable")
for arg, item in it:
if item < current_min:
current_argmin = arg
current_min = item
if return_value:
return current_argmin, current_min
else:
return current_argmin
def argmin(iterable):
return argmin2(enumerate(iterable))
In [52]:
def cuthill_mckee(mat):
"""Return a Cuthill-McKee ordering for the given matrix.
See (for example)
Y. Saad, Iterative Methods for Sparse Linear System,
2nd edition, p. 76.
"""
# this list is called "old_numbers" because it maps a
# "new number to its "old number"
old_numbers = []
visited_nodes = set()
levelset = []
all_nodes = set(xrange(len(mat)))
while len(old_numbers) < len(mat):
if not levelset:
unvisited = list(all_nodes - visited_nodes)
if not unvisited:
break
start_node = unvisited[
argmin(degree(mat, node) for node in unvisited)]
visited_nodes.add(start_node)
old_numbers.append(start_node)
levelset = [start_node]
next_levelset = set()
levelset.sort(key=lambda row: degree(mat, row))
#print levelset
for node in levelset:
row = mat[node]
neighbors, = np.where(row)
for neighbor in neighbors:
if neighbor in visited_nodes:
continue
visited_nodes.add(neighbor)
next_levelset.add(neighbor)
old_numbers.append(neighbor)
levelset = list(next_levelset)
return np.array(old_numbers, dtype=np.intp)
In [53]:
cmk = cuthill_mckee(A)
In [58]:
P = np.eye(len(A))[cmk[::-1]]
A_reordered = np.dot(np.dot(P, A), P.T)
pt.spy(A_reordered, marker=",", precision=prec)
Out[58]:
In [59]:
L = la.cholesky(A_reordered)
print "%d non-zeros" % len(np.where(np.abs(L) > prec)[0])
pt.spy(L, marker=",", precision=prec)
Out[59]:
In []: