numpy中的二进制编码十进制dtype答案

【问题标题】：Binary Coded Decimal dtype in numpynumpy中的二进制编码十进制dtype
【发布时间】：2020-04-18 02:22:47
【问题描述】：

我在 numpy 中使用 matmult（即点积）将一个大型二维矩阵 (dtype=np.bool) 乘以一个一维向量 (dtype=np.uint32)。

np.matmul(matrix, vector, out=vector)

它工作正常，但我在做较大的计算时总是内存不足。

让我印象深刻的一件事是点积的结果向量 - 碰巧 - 我只关心返回的 units - 例如，在返回的结果向量中 - 比如说，一个包含整数 1234，只有 4 个重要，对于 36，只有 6 个重要，依此类推...

这有点远，但与二进制整数溢出时无缝翻转的方式大致相同 - 例如 int8 将按如下方式递增：254、255、256、0、1、2... .

我想知道是否有一种方法可以创建一个只存储一个半字节（如果不支持半字节，则为一个字节）的 dtype，以便 在任何算术运算中只携带十进制单位?

这对于常规二进制编码几乎肯定是不可能的，因为十进制单位存储在 2 的所有幂中。但是，如果 numpy 有一个 BCD 编码的 dtype（或一种有效构造的方法），那么我可能只能存储任何算术运算的 LSB 并且仍然完美地跟踪单元，默默地丢弃每个算术运算上的其他字节；类似于 numpy 的 int 类型溢出的示例。

我知道我可以将 BCD 向后和向前转换为二进制 - 但这没有抓住重点 - 整个计算必须在 BCD 中完成才能工作。任何转换都只需要更多内存。

无论 dtype 被创建来存储 matmult 的结果，必须足够大以存储 length(vector)*max(vector)*max(matrix[row]) - 这通常是一个无符号的32 位数字（特别是对于我的问题，它是 522659*9*1）...对于uint16 来说太大了；但在此之后，我立即使用 (result_vector % 10) 丢弃大部分结果，这将可以存储在 8 位无符号 dtype 中。

内存浪费是相当大的（分析显示使用uint32 作为结果意味着在我的情况下，如果将结果限制为uint8，则需要~1TB 的内存，结果可以存储在~254gb 中）。

那么有什么方法可以通过限制输入的类型来丢弃计算结果 - 使用 BCD 或其他方式？测试表明，如果我将输入向量设为int8，则计算会毫无问题地继续进行，但如上所述会翻转 - 所以它可以使用正确的类型。

但是，我的猜测是我必须完全实现 bcd 类型及其所有操作才能做到这一点？还是实现我自己的自定义矩阵计算？

我很乐意这样做，但想先检查一下我没有错过任何技巧！

最后一件事 - 分析表明 scipy.spare 矩阵无法充分利用矩阵中的零点，因此使用此技巧不会节省内存。索引的成本大于节省的成本，并且比常规 numpy 使用更多的内存。

我已经研究过使用结构化的 dtypes 和视图，这似乎是我正在寻找的东西，但我认为它们都不符合这里的要求。

非常感谢任何想法。

【问题讨论】：

尝试在乘法之前取所有 mod 10。还可以尝试逐行进行点积，一次只占用 O(|row|) 内存。
mod 10 必须发生在矩阵和向量相乘之后。虽然我可以逐行进行（即在每个结果元素上创建它）。
大约是多少。矩阵的形状以及非零条目的百分比有多大？通常二进制矩阵应该占用大部分内存（每个条目 1 个字节，而 uint32 的每个条目 4 个字节）。但是在内存中压缩矩阵也不是太难……
矩阵是三角形的——右上角为1，左下角为0。
所以整个上三角完全被1填满，而下半部分只有0？这将大大简化这个问题。

标签： python python-3.x numpy matrix scipy

【解决方案1】：

与三角布尔矩阵的乘法

在这种情况下，根本不需要巨大的布尔矩阵（522659,522659 需要 273 GB）。 273GB RAM 用于一个简单的算法而不使用内存压缩或其他技巧，2MB RAM 就足够了。

示例

import numpy as np
import numba as nb

@nb.njit(fastmath=True,parallel=True)
def mult_tri(vector):
    out=np.empty(vector.shape[0],dtype=np.uint8)
    for i in nb.prange(vector.shape[0]):
        #accumulator is uint64 to prevent overflow
        acc=np.uint64(0)
        for j in range(i,vector.shape[0]):
            acc+=vector[j]
        out[i]=acc%10
    return out

时间安排

vector = np.random.randint(0,9,522659).astype(np.uint8)
%timeit res_1=mult_tri(vector)
#17.4 s ± 234 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

【讨论】：

【解决方案2】：

您可以使用uint8 类型和一些模除法来减少内存：

import numpy as np

matrix = np.random.choice([True,False],size=(10000,10000))
vector = np.random.randint(0,10000,10000).astype('uint32')

def func1(matrix,vector):
    z = np.empty(1,dtype='uint8')
    v = np.empty(vector.shape[0],dtype='uint8')
    for i,row in enumerate(matrix):
        z = np.tensordot(row,vector,axes=(-1,-1))
        v[i] = z%10
    return v

def func2(matrix,vector):
    z = np.empty(1,dtype='uint8')
    v = np.empty(vector.shape[0],dtype='uint8')
    for i,row in enumerate(matrix):
        np.matmul(row,vector,out=z)
        v[i] = z%10
    return v

matmul 在这种情况下工作得更快。时机明智：

%timeit func1(matrix,vector)
668 ms ± 3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit func2(matrix,vector)
418 ms ± 1.18 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

尝试计算每个总和以迭代地处理uint8 类型的对象是可行的，但是大多数情况下我认为这不值得——这里的中间对象可能是大整数，但对于func2 函数，一次只生成一个一个。因此，即使一行的 matmul 和向量是一个很大的数，一次也只有一个存储在内存中。

numpy 实际上可能会在这里迭代地写入z，在这种情况下，问题没有实际意义——如果有问题，可能值得查看源代码来仔细检查。

【讨论】：