I'm the author. This is for readers interested in spectral methods, structured linear algebra, or interpretable tabular models.
The post looks at restricting learned symmetric matrix pencils, which serve as a nice intermediate between fully interpretable linear models and fully opaque neural networks, to tridiagonal form. That makes the eigensolve much cheaper, still leaves enough expressiveness to fit nontrivial functions, and lets me wire SciPy's `eigh_tridiagonal` into PyTorch autograd for training.
I also compare against the earlier dense version on California Housing: in the 45x45 example, the tridiagonal run takes about 3-4 minutes on CPU versus about 31 minutes for the earlier dense run on an L4 GPU, with similar error in that experiment. The post includes code, plots, and a Colab notebook.
Having developed in C++ a large portion of my career up until a 5+ years ago, it's nice to see C++ went a long way.
However one thing bothers me - error messages. Have the compilers evolved? Eventually, all this stuff translates to templates under the hood, and I remember always getting the compiler error far away from the place I've actually made the mistake.
I am the author of that post series... I see it did some noise :)
I partially agree with your claims. It indeed is an issue if you take high powers of numbers outside of the [-1, 1] interval, but I completely disagree with your assertion that this is the main issue.
Take the famous LogSumExp function as an example - it's used everywhere in machine learning. Nobody says it's a "big no no", even though a naive implementation blows up because you exponentiate large numbers. And if you want to back-prop through it - well, it's even worse. But there is also a "right" way to use it:
LogSumExp(x_1, ..., x_n) = LogSumExp(x_1 - m, ..., x_n - m) + m
where m=max(x_1, ..., x_n)
But the point is that you don't even care! We have good implementations in library functions, such as toch.logsumexp, or scipy.special.logsumexp, and we don't need to worry about it. We never have to compute LogSumExp ourselves and get into the numerical issues.
With polynomials it's the same. Take the Legendre basis. Over an arbitrary interval [a, b], you have the formula
P_n(x) = (2n+1)(x - (a + b)) P_n(x) / (b - a) + n P_n(x)
It's numerically stable. But you don't even care what the formula is - you just call the appropriate function from the np.polynomial.legendre package.
The difference between polynomials and LogSumExp stems from education - people are taught that LogSumExp is "normal" but polynomials are a "big no no". Most people know that there is a library function for LogSumExp, but don't know there are similarly good library functions for the "right" polynomial bases.
And why I believe that the basis is the main issue? Well, it's because even if you use [-1, 1], the standard basis doesn't work with high degrees - it has an extremely badly conditioned vandermonde matrix, and regularization can't help you get a reasonable polynomial. I show it in several posts in the series. The problem with large powers of small numbers is as bad as with large numbers. And the issue with the standard basis is that you have no choice, but to compute these high powers - that's how the basis is defined.
Having said that, I will obviously take your input into account, and say something about high powers in the post. Explain why I believe it's exactly like LogSumExp, and why I believe it's the basis (rather than floating point errors) that is the issue. You are correct that I wrongly assumed it's obivious, and it is not.
Like with any object, oftentimes working with it by definition is very different with how you work with it computationally. But I believe that once you have good library function, you don't really care.
The post looks at restricting learned symmetric matrix pencils, which serve as a nice intermediate between fully interpretable linear models and fully opaque neural networks, to tridiagonal form. That makes the eigensolve much cheaper, still leaves enough expressiveness to fit nontrivial functions, and lets me wire SciPy's `eigh_tridiagonal` into PyTorch autograd for training.
I also compare against the earlier dense version on California Housing: in the 45x45 example, the tridiagonal run takes about 3-4 minutes on CPU versus about 31 minutes for the earlier dense run on an L4 GPU, with similar error in that experiment. The post includes code, plots, and a Colab notebook.