Преобразование типов в numpy, защита от переполнения

одно из решений - np.linalg.norm это - норма вектора, корень из суммы квадратов, как раз то что надо. Внутри оно использует float, и выглядит удовлетворительным

samples = np.frombuffer(indata, dtype=np.int16)[::scale] #int16
l=len(samples)
en = np.linalg.norm(samples)/np.sqrt(l) #float
rms = np.int16(en) #int16

numpy 2.3.3

Если функция не запрашивает тип результата, но работает правильно, то скорее всего где-то внутри она создаёт копии исходных данных, преобразованных к приемлемому типу.

Например, функция numpy.linalg.norm с первых строк проверяет тип данных, и если это целые числа, то они будут преобразованы в numpy.float64. То есть, применяя linalg.norm к целочисленному вектору, мы всё равно расходуем дополнительную память и время на преобразование типов, просто от нас это скрыто:

    # выдержка из кода numpy.linalg.norm, где x - переданный вектор
    ...
    if not issubclass(x.dtype.type, (inexact, object_)):
        x = x.astype(float)
    ...

Но замечу, что для вектора типа numpy.float16 скрытого преобразования в 8-байтовый float не произойдёт, и мы рискуем получить бесконечность из-за переполнения. То есть, для float16 нам нужно самостоятельно преобразовать тип, чтобы избежать ошибки.

Практические выводы я бы сделал такие:

• Если вычислений много, то лучше один раз самому сделать приведение к приемлемому типу, чтобы избежать многочисленного скрытого преобразования.
• Если преобразование нужно только в одном месте, то можно попробовать делегировать как можно больше работы внутрь скомпилированного кода NumPy.

Что касается второго пункта, то я бы обратил внимание на подходящие функции, которые принимают параметр dtype=..., см. Universal functions:

dtype

Overrides the DType of the output arrays the same way as the signature. This should ensure a matching precision of the calculation. The exact calculation DTypes chosen may depend on the ufunc and the inputs may be cast to this DType to perform the calculation.

---

Переопределяет DType выходных массивов так же, как параметр signature. Это должно обеспечить соответствующую точность вычислений. Выбранные точные типы вычислений DTypes могут зависеть от ufunc, и входные данные могут быть преобразованы в этот DType для выполнения вычислений.

Доступные комбинации типов входных и выходных данных можно посмотреть с помощью свойства types. Например:

>>> print(np.vecdot.types)

['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I',
 'll->l', 'LL->L', 'qq->q', 'QQ->Q', 'ee->e', 'ff->f', 'dd->d',
 'gg->g', 'FF->F', 'DD->D', 'GG->G', 'OO->O']

Здесь две буквы слева - код типа входных данных, буква справа - выходных, код типов можно посмотреть в документации или через свойтво char, например numpy.dtype(numpy.float64).char равен 'd'.

Внутри такая логика: по заданному типу результата делается приведение входных данных к подходящему типу, после чего вызывается работающая с этим типом версия функции. Продолжая пример с vecdot, для результата float применяется правило 'dd->d', т.е. входные данные тоже преобразуются к типу float перед вызовом нужной версии основной функции.

Из функций, которые можно применить для решения исходной задачи, я бы выделил multiply, matmul и vecdot. Нюанс в том, что они принимают два аргумента. Моих знаний Си не хватает, чтобы сказать, будут ли оба аргумента приведены к заданному типу независимо друг от друга, когда мы передаём один и тот же массив дважды. Я бы исходил из того, что в работе с одним массивом будет лучше явно преобразовать его тип и подставить результат дважды как аргумент. Например:

np.sqrt(np.vecdot((_:=sample.astype(float)), _) / len(sample))

вместо

np.sqrt(np.vecdot(sample, sample, dtype=float) / len(sample))

Что касается скорости работы, то на моей стороне результат на разных функциях соизмерим. Но если среди рассмотренных вариантов что-то выделить, то я бы обратил внимание на vecdot и matmul:

import numpy as np
from timeit import Timer
from functools import partial
import matplotlib.pyplot as plt
import seaborn as sns

code = {
    'norm':           'np.linalg.norm(sample.astype(float))/np.sqrt(len(sample))',
    'vecdot _':       'np.sqrt(np.vecdot((_:=sample.astype(float)), _) / len(sample))',
    'vecdot dtype':   'np.sqrt(np.vecdot(sample, sample, dtype=np.float64) / len(sample))',
    'multiply _':     'np.sqrt(np.multiply((_:=sample.astype(float)), _).sum() / len(sample))',
    'multiply dtype': 'np.sqrt(np.multiply(sample, sample, dtype=np.float64).sum() / len(sample))',
    'matmul _':       'np.sqrt(np.matmul((_:=sample.astype(float)), _) / len(sample))',
    'matmul dtype':   'np.sqrt(np.matmul(sample, sample, dtype=np.float64) / len(sample))',
}

# Sanity check
test_globals={'np': np, 'sample': np.random.randint(-32768, 32768, size=16_000, dtype=np.int16)}
test_eval = partial(eval, globals=test_globals)
assert np.allclose((_:=[*map(test_eval, code.values())])[:-1], _[1:])
assert eval(code['norm'], globals={'np': np, 'sample': np.array([1, 1, 1])}) == 1

# timeit
r, n  = 50, 100
setup = 'sample = np.random.randint(-32768, 32768, size=16_000, dtype=np.int16)'
stats = {name: Timer(stmt, setup, globals={'np': np}).repeat(r, n) for name, stmt in code.items()}

# plot
fig, ax = plt.subplots(figsize=(15, 5))
ax.tick_params(axis='both', which='both', length=0)
custom_colors = ['OrangeRed', 'Orange', 'DarkOrange', 'Lime', 'LimeGreen', 'DeepSkyBlue', 'DodgerBlue',]
sns.violinplot(stats, orient='h', palette=custom_colors, ax=ax) # showfliers=False