That's an interesting question. I don't have too much experience, but here's my two cents.
For matrix function approximations, loss of orthogonality matters less than for eigenvalue computations. The three-term recurrence maintains local orthogonality reasonably well for moderate iteration counts. My experiments [1] show orthogonality loss stays below $10^{-13}$ up to k=1000 for well-conditioned problems, and only becomes significant (jumping to $10^{-6}$ and higher) around k=700-800 for ill-conditioned spectra. Since you're evaluating $f(T_k)$ rather than extracting individual eigenpairs, you care about convergence of $\|f(A)b - x_k\|$, not spectral accuracy. If you need eigenvectors themselves or plan to run thousands of iterations, you need the full basis, and the two-pass method won't help. Maybe methods like [2] would be more suitable?
Thanks! I used perf to look at cache miss rates and memory bandwidth during runs. The measurements showed the pattern I expected, but I didn't do a rigorous profiling study (different cache sizes, controlled benchmarks across architectures, or proper statistical analysis).
This was for a university exam, and I ran out of time to do it properly. The cache argument makes intuitive sense (three vectors cycling vs. scanning a growing n×k matrix), and the timing data supports it, but I'd want to instrument it more carefully in the future :)
The blog post is a simplification of the actual work; you can check out the full report here [1], where I also reference the literature about this algorithm.
On the cache effects: I haven't seen this "engineering" argument made explicitly in the literature either. There are other approaches to the basis storage problem, like the compression technique in [2]. Funny enough, the authors gave a seminar at my university literally this afternoon about exactly that.
I'm also unfamiliar with randomised algorithms for numerical linear algebra beyond the basics. I'll dig into that, thanks!
On the BLAS point, let me clarify what I meant by "wall": when you call BLAS from Rust, you're essentially making a black-box call to pre-compiled Fortran or C code. The compiler loses visibility into what happens across that boundary. You can't inline, can't specialise for your specific matrix shapes or use patterns, can't let the compiler reason about memory layout across the whole computation. You get the performance of BLAS, sure, but you lose the ability to optimise the full pipeline.
Also, Rust's compilation model flattens everything into one optimisation unit: your code, dependencies, all compiled together from source. The compiler sees the full call graph and can inline, specialise generics, and vectorise across what would be library boundaries in C/C++. The borrow checker also proves at compile time that operations like our pointer swaps are safe and that no aliasing occurs, which enables more aggressive optimisations; the compiler can reorder operations and keep values in registers because it has proof about memory access patterns. With BLAS, you're calling into opaque binaries where none of this analysis is possible.
My point is that if the core computation just calls out to pre-compiled C or Fortran, you lose much of what makes Rust interesting for numerical work in the first place. That's why I hope to see more efforts directed towards expanding the Rust ecosystem in this area in the future :)
Hi there, thanks! I started doing this for a university exam and got carried away a bit.
Regarding Rust for numerical linear algebra, I kinda agree with you. I think that theoretically, its a great language for writing low-level "high-performance mathematics." That's why I chose it in the first place.
The real wall is that the past four decades of research in this area have primarily been conducted in C and Fortran, making it challenging for other languages to catch up without relying heavily on BLAS/LAPACK and similar libraries.
I'm starting to notice that more people are trying to move to Rust for this stuff, so it's worth keeping an eye open on libraries like the one that I used, faer.
For many applications, casting to u16 and wasting 6 bits is perfectly fine. The "trouble" is only worth it when you're operating at a scale where those wasted bits add up to gigabytes.
This is common in fields like bioinformatics, search engine indexing, or implementing other succinct data structures. In these areas, the entire game is about squeezing massive datasets into RAM to avoid slow disk I/O. Wasting 6 out of 16 bits means your memory usage is almost 40% higher than it needs to be. That can be the difference between a server needing 64GB of RAM versus 100GB.
On top of that, as I mentioned in another comment, packing the data more tightly often makes random access faster than a standard Vec, not slower. Better cache locality means the CPU spends less time waiting for data from main memory, and that performance gain often outweighs the tiny cost of the bit-fiddling instructions.
It's not about poorer instructions; a get_unchecked on a Vec<u8> is just a single memory access, which is as good as it gets. The difference is likely down to cache locality effects created by the benchmark loop itself.
The benchmark does a million random reads. For the FixedVec implementation with bit_width=8, the underlying storage is a Vec<u64>. This means the data is 8x more compact than the Vec<u8> baseline for the same number of elements.
When the random access indices are generated, they are more likely to fall within the same cache lines for the Vec<u64>-backed structure than for the Vec<u8>. Even though Vec<u8> has a simpler access instruction, it suffers more cache misses across the entire benchmark run.
The other replies nailed it, but I'll add my two cents.
Vec<u8> may be the right call most of the time for most use cases. This library, however, is for when even 8 bits compared to 4 is too much. Another example, if all your values fit in 9 bits, you'd be forced to use Vec<u16> and waste 7 bits for every single number. With this structure, you just use 9 bits. That's almost a 2x space saving. At scale, that makes a difference.
For floats, you'd use fixed-point math, just like the other replies said. Your example of +/- 10.0 with 1e-6 precision would mean multiplying by 1,000,000 and storing the result as a 25-bit signed integer. When you read it back, you just divide. It's a decent saving over a 32-bit float.
That's a good catch!! Thank you, you are right. I (incorrectly) assumed that a single u64 could capture the entire bit_width-value read starting from byte_pos. However, as you said, this assumption breaks for some large bit widths.
BEXTR basically does the same thing, yes. I'm sticking with the portable shift-and-mask, though. My bet is that LLVM is smart enough to see that pattern and emit a BEXTR on its own when the target supports BMI1.
Using the intrinsic directly would also kill portability. I'd need #[cfg] for ARM and runtime checks for older x86 CPUs, which adds complexity for a tiny, if any, gain. The current code just lets the compiler pick the best instruction on any architecture.
For matrix function approximations, loss of orthogonality matters less than for eigenvalue computations. The three-term recurrence maintains local orthogonality reasonably well for moderate iteration counts. My experiments [1] show orthogonality loss stays below $10^{-13}$ up to k=1000 for well-conditioned problems, and only becomes significant (jumping to $10^{-6}$ and higher) around k=700-800 for ill-conditioned spectra. Since you're evaluating $f(T_k)$ rather than extracting individual eigenpairs, you care about convergence of $\|f(A)b - x_k\|$, not spectral accuracy. If you need eigenvectors themselves or plan to run thousands of iterations, you need the full basis, and the two-pass method won't help. Maybe methods like [2] would be more suitable?
[1] https://github.com/lukefleed/two-pass-lanczos/raw/master/tex...
[2] https://arxiv.org/abs/2403.04390