Computing the kernel of a tridiagonal matrix. Suppose $M = M(\lambda)$ is a symmetric  $n\times n$ tridiagoal matrix that depends polynomially on a real number $\lambda$ in the sense that
$$
M(\lambda) = A_0 + \lambda A_1 + \dots + \lambda^4 A_4, 
$$
where $A_i$'s are known tridiagonal matrices.

What is a good numerical algorithm that would allow us to compute $\ker M(\lambda)$ ?

An obvious method would be to first solve the equation
$$
\det M(\lambda) = 0
$$
for all possible $\lambda\in\Bbb R$ then find the corresponding $v_\lambda \in \ker M(\lambda)$. I don't know if this method is efficient or stable enough, or what would be the best way to solve the equation $\det M(\lambda) = 0$. I would be very thankful if anyone knows an efficient algorithm to solve this equation, or skip it to compute $\ker M(\lambda)$ directly.
P.S. I was told that this problem arose from a certain discretization of some eigenvalue problem related to a PDE.
 A: I assume that $M(\lambda)$ is not singular unless for some $\lambda = \lambda_i$ where $\det M(\lambda_i) = 0$ and $\dim \ker M(\lambda_i) = 1$. I also assume that all roots of $\det M(\lambda) = 0$ are simple.
Solving $\det M(\lambda) = 0$ directly is hard, since $\det M(\lambda)$ would a polynomial of $4n$ degree where $n$ is the number of rows in $M$.
The following numerical approach may be used to find the $\lambda$ and the kernel simultaneously:

*

*Choose some segment $[a, b]$. The method will search for eigenvalues within this segment.

*Define function $\sigma(\lambda)$ which is the signature of $M(\lambda)$. The signature basically is the number of positive eigenvalues minus the number of negative eigenvalues. Sylvester's law of inertia states that you don't really need to compute the eigenvalues, but instead you may compute the leading minors of the matrix. And those minors are easy to compute if you have factorized your matrix using LU or LDL (preffered) decomposition.

*$\sigma(\lambda)$ is a piecewise-constant function jumping by 2 every time when $\lambda$ crosses $\lambda_i$. This is because when $\lambda = \lambda_i$ some eigenvalue of $M(\lambda)$ changes its sign.

*Use bisection method to bracket each jump location (which are your sought values $\lambda_i$)

Note that compute the kernel is rather simple if you already have the decomposition.
$$
M(\lambda_i) = LDL^\top\\
L = \begin{pmatrix}
1 \\
l_{2,1} & 1\\
\vdots & \ddots & 1\\
l_{r,1} & \cdots & l_{r,r-1} & 1\\
\ast & & \cdots & & \ast\\
\vdots & & & & & \ddots\\
\ast & & & \cdots & & & \ast
\end{pmatrix}\\
D = \begin{pmatrix}
d_1 \\
& d_2\\
&& \ddots\\
&&& d_r \\
&&&& \ast \\
&&&&& \ddots\\
&&&&&& \ast
\end{pmatrix}
$$
Here is assumed that decomposition breaks at row $r$ with $|d_r| < \epsilon$. Not computed values marked as $\ast$.
Solving
$$
LDL^\top x = 0
$$
reduces to
$$
\begin{pmatrix}
1 & l_{2,1} & \cdots & l_{r,1} & \ast & \cdots & \ast\\
 & 1 & \ddots & \vdots &  & & \\
&& 1 & l_{r,r-1} & \vdots & & \\
&&& 1 &  & & \vdots\\
&&&& \ast\\
&&&&& \ddots\\
&&&&&& \ast\\
\end{pmatrix}
\begin{pmatrix}
x_1\\
\vdots\\
x_r\\
x_{r+1}\\
\vdots\\
x_n
\end{pmatrix} = 
\begin{pmatrix}
0\\
\vdots\\
1\\
0\\
\vdots\\
0
\end{pmatrix}
=
\mathbf e_r
$$
Plugging $x_r = 1, x_{r+1} = \cdots = x_n = 0$ and solving the system via back substitution gives the desired kernel vector. Note that for tridiagonal case the $L$ matrix would be bidiagonal and the factorization and solution steps have $O(n)$ complexity.
Here is a Jupyter notebook with a proof-of-concept implementation.
