I am trying to efficiently solve the following linear system in Python 3.6:
b = T x
where T is an N x N Toeplitz matrix and x, b are N x N matrices. It is much slower than both scipy.linalg.solve and numpy.linalg.solve (see the output from my test script for N=500).
A possible explanation is that scipy.linalg.solve_toeplitz uses a for-loop for matrix inputs to solve the individual linear systems with x, b being vectors (solving each column at a time). Unfortunately, i was unable to make the benchmark work from the original pull request adding this function.
Reproducing code example:
"""This script serves to benchmark several methods for solving a
linear equation involving a Toeplitz matrix and determine which one is
faster.
We consider the equation: b = T x, where T is Toeplitz, b and x are matrices and we wish to
solve for x.
"""
import numpy as np
import numpy.linalg
import scipy.linalg
import time
N = 500
np.random.seed(1)
# Construct random vectors/matrices
x = np.random.rand(N, N)
c,r = np.random.rand(2, N)
T = scipy.linalg.toeplitz(c, r)
b = T.dot(x)
# Validate various solutions
x_sol = []
x_sol.append(('numpy.linalg.solve(T, b)', numpy.linalg.solve(T, b)))
x_sol.append(('scipy.linalg.solve(T, b)', scipy.linalg.solve(T, b)))
x_sol.append(('scipy.linalg.solve_toeplitz((c, r), b)', scipy.linalg.solve_toeplitz((c, r), b)))
for solution in x_sol:
print("Method: {} - error: {}".format(solution[0], numpy.linalg.norm(solution[1] - x)))
# Time the solutions
x_time = []
for solution in x_sol:
t_start = time.time()
eval(solution[0])
t_end = time.time()
x_time.append(t_end - t_start)
print("Timings:")
print(x_time)
Test output:
Method: numpy.linalg.solve(T, b) - error: 5.996096106604476e-11
Method: scipy.linalg.solve(T, b) - error: 5.5520852583734864e-11
Method: scipy.linalg.solve_toeplitz((c, r), b) - error: 2.761917504645538e-08
Timings:
[0.008001089096069336, 0.44705653190612793, 0.2655336856842041]
Scipy/Numpy/Python version information:
'0.19.1', '1.13.3', sys.version_info(major=3, minor=6, micro=3, releaselevel='final', serial=0)
I am trying to efficiently solve the following linear system in Python 3.6:
b = T x
where T is an N x N Toeplitz matrix and x, b are N x N matrices. It is much slower than both scipy.linalg.solve and numpy.linalg.solve (see the output from my test script for N=500).
A possible explanation is that scipy.linalg.solve_toeplitz uses a for-loop for matrix inputs to solve the individual linear systems with x, b being vectors (solving each column at a time). Unfortunately, i was unable to make the benchmark work from the original pull request adding this function.
Reproducing code example:
Test output:
Method: numpy.linalg.solve(T, b) - error: 5.996096106604476e-11
Method: scipy.linalg.solve(T, b) - error: 5.5520852583734864e-11
Method: scipy.linalg.solve_toeplitz((c, r), b) - error: 2.761917504645538e-08
Timings:
[0.008001089096069336, 0.44705653190612793, 0.2655336856842041]
Scipy/Numpy/Python version information: