I'm trying to optimize a piece of code that solves a large sparse nonlinear system using an interior point method. During the update step, this involves computing the Hessian matrix `H`

, the gradient `g`

, then solving for `d`

in `H * d = -g`

to get the new search direction.

The Hessian matrix has a symmetric tridiagonal structure of the form:

A.T * diag(b) * A + C

I've run `line_profiler`

on the particular function in question:

```
Line # Hits Time Per Hit % Time Line Contents
==================================================
386 def _direction(n, res, M, Hsig, scale_var, grad_lnprior, z, fac):
387
388 # gradient
389 44 1241715 28220.8 3.7 g = 2 * scale_var * res - grad_lnprior + z * np.dot(M.T, 1. / n)
390
391 # hessian
392 44 3103117 70525.4 9.3 N = sparse.diags(1. / n ** 2, 0, format=FMT, dtype=DTYPE)
393 44 18814307 427597.9 56.2 H = - Hsig - z * np.dot(M.T, np.dot(N, M)) # slow!
394
395 # update direction
396 44 10329556 234762.6 30.8 d, fac = my_solver(H, -g, fac)
397
398 44 111 2.5 0.0 return d, fac
```

Looking at the output it's clear that constructing `H`

is by far the most costly step - it takes considerably longer than actually solving for the new direction.

`Hsig`

and `M`

are both CSC sparse matrices, `n`

is a dense vector and `z`

is a scalar. The solver I'm using requires `H`

to be either a CSC or CSR sparse matrix.

Here's a function that produces some toy data with the same formats, dimensions and sparseness as my real matrices:

```
import numpy as np
from scipy import sparse
def make_toy_data(nt=200000, nc=10):
d0 = np.random.randn(nc * (nt - 1))
d1 = np.random.randn(nc * (nt - 1))
M = sparse.diags((d0, d1), (0, nc), shape=(nc * (nt - 1), nc * nt),
format='csc', dtype=np.float64)
d0 = np.random.randn(nc * nt)
Hsig = sparse.diags(d0, 0, shape=(nc * nt, nc * nt), format='csc',
dtype=np.float64)
n = np.random.randn(nc * (nt - 1))
z = np.random.randn()
return Hsig, M, n, z
```

And here's my original approach for constructing `H`

:

```
def original(Hsig, M, n, z):
N = sparse.diags(1. / n ** 2, 0, format='csc')
H = - Hsig - z * np.dot(M.T, np.dot(N, M)) # slow!
return H
```

Timing:

```
%timeit original(Hsig, M, n, z)
# 1 loops, best of 3: 483 ms per loop
```

Is there a faster way to construct this matrix?