Gilbert Strang explains the matrix square root of the second difference matrix here.
In particular,
import numpy
import scipy.linalg
def f(p):
return (4 / numpy.pi) / ((1 - 2*p) * (1 + 2*p))
N = 10000
K_tri = numpy.eye(N) - numpy.eye(N, k=1)
K = K_tri + K_tri.T
s = numpy.arange(N)
T = scipy.linalg.toeplitz(f(s))
H = -scipy.linalg.hankel(f(s+2), f(N+1-s))
K_sqrt = T + H
print(numpy.max(numpy.abs(numpy.dot(K_sqrt, K_sqrt) - K)))
7.90678633678e-12