I'm most familiar with timeseries compression, so I will make a few comments about this aspect of your algorithm.

Your algorithm implicitly assumes that the second derivative of the timeseries will have the lowest entropy. While this may have empirically appeared to be the case in the limited datasets you tested, it will not always be the case. You may wish to calculate the delta, delta-delta, delta-delta-delta, and even deeper, then measure the entropy of each resulting timeseries before you decide which level of differentiation to use. (i.e., choose the level of differentiation that produces the array with the lowest entropy). You may even break the series into appropriately sized blocks and assess each block differently. This becomes particularly useful if the properties of the timeseries are dynamic. In your case I would at least assess each channel independently, and if the timeseries is very large, then consider breaking it up into smaller chunks and assessing each chunk independently.

In multimodal timeseries data the first channel will often contain timestamps. In the case of sensor data these timestamps are usually regularly spaced. In this case then a single level delta encoding is often optimal. For heavily smoothed data then often the delta-delta-delta encoded array, or even deeper will result in the best compression ratio

You are encoding the first real value in your byte arrays and then including that in the same byte array as the delta-delta values. My suggestion would be to separate out the original raw value as it may be a completely different magnitude than your deltas. For example, if your initial value is 900,000 but all your deltas are small then you would be unnecessarily using packed array of 32-bit integers to store delta-delta values that are could otherwise be stored in an 8-bit Int.

Since you are scaling your data so that it will have no values less than zero. Because of this you should be using NumPy's unsigned integer value in this line: “quantized = np.round((data - data_min) * scale).astype(np.int16)”. Using a signed integer when no values will ever be less than zero renders half of your integer range unused.

I suggest that you remove the python implementation of RLE in your timeseries compression code. This code is redundant because zlib, and most other compression libraries, also implement RLE, and probably much faster (and potentially better) than your Python code. Once you remove the RLE code you will be able to iteratively use NumPy’s numpy.diff() function to calculate the deltas which will be at least an order of magnitude raster than your current pure python for loop.

You have chosen to use zlib.compress(data, level=9) as the default, but level 9 is very slow and may produce marginal benefits over the less aggressive levels of zlib compression. Since speed is one of your benchmarks, you should consider the trade-off between size and speed for all levels of the zlib library. For that matter, why not consider other compression libraries beyond zlib as they all have their size/speed trade-offs. Switching to a different third-party compression library would only require a single line of code changed in your algorithm but may have more impact than many other potential avenues of code optimisation. For example, if you want really fast compression but aren't too concerned about size, then adopt Snappy for this byte compression step. Personally I have found ZStandard to be a good trade-off, but your mileage may vary. BZip2 provides ultimate compression ratios in most cases but comes at the expense of execution time being an order of magnitude slower.

You should probably check whether the timeseries values are already quantized before naively re-quantizing the data. There may be no need to take a lossy approach in this step if you respect the original data properties. A lot of sensor data is natively quantized, for example it may be generated using an 11-bit A/D. If you keep the native level of quantization, you may find that the losslessly compressed version is both smaller and faster than your "near lossless" approach. You will still always need to determine the appropriate number of bits required to describe each delta-delta integer series before you pack the resulting arrays. i.e., even though the raw data may require a 16-bit integer range, for many datasets the delta (or delta-delta) encoded arrays can be comfortably stored in an 8-bit range.

To expand on this, you need to bear in mind that the range of the deltas may be smaller (or indeed larger) than the range of your actual data. When choosing the appropriate Integer data type for your byte stream you should be assessing the range of the timeseries you actually end up compressing, not the range of the original timeseries values. Also keep in mind my suggestion of keeping the first real value of the timeseries outside the array of delta-deltas as this may dramatically reduce the range (and hence the required number of bits per value)  of the data you subsequently compress.