标签: linear-algebra

我试图找出在python中找到稀疏对称和实矩阵行列式的最快方法.使用scipy sparse模块,但真的很惊讶没有行列式功能.我知道我可以使用LU分解来计算行列式,但是没有看到一个简单的方法来执行它,因为返回scipy.sparse.linalg.splu是一个对象并且实例化一个密集的L和U矩阵是不值得的 - 我不妨sp.linalg.det(A.todense())在哪里做A我的scipy稀疏矩阵.

我也有点惊讶为什么其他人没有面对scipy中有效决定因素计算的问题.如何使用splu计算行列式？

我看着pySparse和scikits.sparse.chlmod.后者对我来说现在不实用 - 需要软件包安装,也不确定在我遇到麻烦之前代码的速度有多快.有解决方案吗 提前致谢.

我有一个(nxn)对称矩阵A和一个(nx 1)向量B.基本上,我只需要解决Ax = b for x.问题是A可能会很大.所以我正在寻找最有效的C++求解线性方程的算法.我查看了Eigen库.显然它有一个SVD方法,但我被告知它很慢.解决x =逆(A)*b似乎也不是最理想的.uBLAS更快吗？有没有更有效的方法？谢谢.

编辑:矩阵A是正定的而不是稀疏的.

(出于测试目的)我编写了一个简单的方法来计算nxn矩阵的转置

void transpose(const size_t _n, double* _A) {
    for(uint i=0; i < _n; ++i) {
        for(uint j=i+1; j < _n; ++j) {
            double tmp  = _A[i*_n+j];
            _A[i*_n+j] = _A[j*_n+i];
            _A[j*_n+i] = tmp;
        }
    }
}

当使用优化级别O3或Ofast时,我期望编译器展开一些循环,这将导致更高的性能,尤其是当矩阵大小是2的倍数(即,每次迭代可以执行双循环体)或类似时.相反,我测量的恰恰相反.2的权力实际上表明执行时间显着增加.

这些尖峰也是64的固定间隔,间隔128更明显,依此类推.每个尖峰延伸到相邻的矩阵大小,如下表所示

size n  time(us)
1020    2649
1021    2815
1022    3100
1023    5428
1024    15791
1025    6778
1026    3106
1027    2847
1028    2660
1029    3038
1030    2613

我使用gcc版本4.8.2编译但是同样的事情发生在clang 3.5上,所以这可能是一些通用的东西？

所以我的问题基本上是:为什么执行时间周期性增加？是否有一些通用的东西与任何优化选项一起出现(就像clang和gcc一样)？如果是这样的优化选项导致了这个？

这怎么可能如此重要,即使O0版本的程序在512的倍数时优于03版本？

编辑:注意此(对数)图中峰值的大小.转换具有优化的1024x1024矩阵实际上花费的时间与在没有优化的情况下转置1300x1300矩阵一样多.如果这是一个缓存故障/页面错误问题,那么有人需要向我解释为什么内存布局对于程序的优化版本是如此显着不同,它失败的功能为2,只是为了恢复高性能稍大的矩阵.缓存故障是否应该创建更多类似步骤的模式？为什么执行时间会再次下降？(为什么优化会创建以前不存在的缓存错误？)

编辑:以下应该是gcc生成的汇编代码

没有优化(O0):

_Z9transposemRPd:
.LFB0:
    .cfi_startproc
    push    rbp …

在我学会了如何使用之后einsum,我现在正试图了解它是如何np.tensordot工作的.

但是,我有点迷失,特别是关于参数的各种可能性axes.

要理解它,因为我从未练习过张量微积分,我使用以下示例:

A = np.random.randint(2, size=(2, 3, 5))
B = np.random.randint(2, size=(3, 2, 4))

在这种情况下,有什么不同可能np.tensordot,你会如何手动计算？

给定一组带有x,y坐标的几百万个点,快速找到一个位置的前1000个最近点的算法是什么？"快速"在这里意味着家用电脑上大约100毫秒.

蛮力意味着进行数百万次乘法,然后对它们进行排序.虽然一个简单的Python应用程序可以在不到一分钟的时间内完成,但对于交互式应用程序来说仍然太长.

点的边界框将是已知的,因此将空间划分为简单网格是可能的.然而,点的分布有些不均匀,所以我怀疑大多数网格方块都是空的,然后突然其中一些将包含大部分点.

编辑:不必确切,实际上可能非常不准确.如果前1000名实际上只是来自前2000名的一些随机点,那就没什么大不了的了.

编辑:点集很少改变.

问候,谷歌搜索该主题带来,例如,MTL,exmat,LAPACK以及此处.我似乎也记得微软研究院发布了一个,但不能把它放在手上.我寻求实际使用(或开发......)其中之一的人的建议,希望在C++中尽可能地获得Matlab经验(尽可能多).谢谢,罗比

我不需要任何花哨的东西; 基本上,如果它会让我做一些基本的Matrix*Matrix的东西,我会很高兴.

为什么BLAS具有gemm矩阵 - 矩阵乘法gemv功能和矩阵向量乘法的单独函数？矩阵向量乘法不是矩阵乘法的一个特例,其中一个矩阵只有一个行/列吗？

TL; DR如何在没有重建的情况下将ATLAS/MKL链接到现有的Numpy.

我用Numpy用大矩阵计算,我发现它非常慢,因为Numpy只使用1个核心进行计算.经过大量的搜索,我发现我的Numpy没有链接到像ATLAS/MKL这样的优化库.这是我的numpy配置:

>>>import numpy as np
>>>np.__config__.show()
blas_info:
    libraries = ['blas']
    library_dirs = ['/usr/lib']
    language = f77
lapack_info:
    libraries = ['lapack']
    library_dirs = ['/usr/lib']
    language = f77
atlas_threads_info:
    NOT AVAILABLE
blas_opt_info:
    libraries = ['blas']
    library_dirs = ['/usr/lib']
    language = f77
    define_macros = [('NO_ATLAS_INFO', 1)]
atlas_blas_threads_info:
  NOT AVAILABLE
openblas_info:
  NOT AVAILABLE
lapack_opt_info:
    libraries = ['lapack', 'blas']
    library_dirs = ['/usr/lib']
    language = f77
    define_macros = [('NO_ATLAS_INFO', 1)]
atlas_info:
  NOT AVAILABLE
lapack_mkl_info:
  NOT AVAILABLE
blas_mkl_info:
  NOT AVAILABLE
atlas_blas_info:
  NOT AVAILABLE
mkl_info: …

是否有任何算法可以从一组3D数据点返回直线方程？我可以找到很多来源,这些来源将给出2D数据集中的线的等式,但没有3D.

谢谢.