我正在尝试解决一组形式为Ax = 0的方程式.A是已知的6x6矩阵,我使用SVD编写了下面的代码,以获得在某种程度上起作用的向量x.答案大致正确,但不足以对我有用,我怎样才能提高计算的精确度？降低低于1.e-4的eps会导致功能失败.
from numpy.linalg import *
from numpy import *
A = matrix([[0.624010149127497 ,0.020915658603923 ,0.838082638087629 ,62.0778180312547 ,-0.336 ,0],
[0.669649399820597 ,0.344105317421833 ,0.0543868015800246 ,49.0194290212841 ,-0.267 ,0],
[0.473153758252885 ,0.366893577716959 ,0.924972565581684 ,186.071352614705 ,-1 ,0],
[0.0759305208803158 ,0.356365401030535 ,0.126682113674883 ,175.292109352674 ,0 ,-5.201],
[0.91160934274653 ,0.32447818779582 ,0.741382053883291 ,0.11536775372698 ,0 ,-0.034],
[0.480860406786873 ,0.903499596111067 ,0.542581424762866 ,32.782593418975 ,0 ,-1]])
def null(A, eps=1e-3):
u,s,vh = svd(A,full_matrices=1,compute_uv=1)
null_space = compress(s <= eps, vh, axis=0)
return null_space.T
NS = null(A)
print "Null space equals ",NS,"\n"
print dot(A,NS)
我正在进行大量的科学编程,并使用Boost.Units提供了非常好的经验,它提供了数量的编译时尺寸分析(即带有单位的标签数量,从而通过经典物理尺寸分析捕获了许多误差)并使用了Eigen 2代表线性代数.
然而,Eigen没有单位的概念,虽然你可以在Eigen的矩阵中设置标量,但是它期望两个量的乘法产生相同的类型,这对于单位来说显然是不正确的.例如,代码如:
using boost::units::quantity;
namespace si = boost::units::si;
Eigen::Matrix< quantity< si::length >, 2, 1 > meter_vector;
quantity< si::area > norm = meter_vector.squaredNorm();
不起作用,即使它在逻辑上是正确的.
有没有支持单位的矩阵库？我知道这在过去很难实现,而且C++ 11 decltype会更容易实现,但是C++ 03和模板专业化肯定是可能的.
我有一个关于使用CUDA对数百个小矩阵进行特征分解的问题.
我需要同时计算数百(例如500)小(64乘64)实对称矩阵的特征值和特征向量.我试图通过Jacobi方法使用国际象棋锦标赛订购来实现它(有关更多信息,请参阅本文(PDF)).
在该算法中,在每个块中定义了32个线程,而每个块处理一个小矩阵,并且32个线程一起工作以使32个非对角线元素膨胀直到收敛.但是,我对它的表现并不十分满意.
我想知道我的问题哪里有更好的算法,即许多64乘64实对称矩阵的特征分解.我想家庭主人的方法可能是更好的选择,但不确定它是否可以在CUDA中有效实施.网上没有很多有用的信息,因为大多数其他程序员更感兴趣的是使用CUDA/OpenCL来分解一个大矩阵而不是很多小矩阵.
我认为这可能是一个简单的数学问题,但我不知道现在发生了什么.
我在网络摄像头上捕捉"标记"的位置,我有一个标记及其坐标列表.四个标记是工作表面的外角,第五个(绿色)标记是小部件.像这样:
替代文字http://i37.tinypic.com/308cjtv.jpg
这是一些示例数据:
我想以某种方式将网络摄像头的小部件位置转换为坐标以显示在屏幕上,其中左上角是0,0而不是98,86并且以某种方式考虑了网络摄像头捕获的扭曲角度.
我甚至会从哪里开始？任何帮助赞赏
我正在尝试将一些Python代码移植到Scala.它大量使用Numpy和Scipy.虽然我发现了许多密集矩阵/线性代数库,它们可以作为NumPy的一个适当的(但不是极好的)替代品,但我还没有找到任何提供我在SciPy中使用的功能的东西.特别是,我正在寻找支持稀疏部分特征分解的库(比如SciPy包装的arpack),然后是SciPy提供的一些简单事物的库(例如直方图).
在Eigen版本中,我使用"真正的"固定大小矩阵和向量,更好的算法(LDLT与uBlas的LU),它在内部使用SIMD指令.那么,为什么它在下面的例子中比uBlas慢？
我确信,我做错了 - Eigen 必须更快,或至少可比.
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/ublas/vector.hpp>
#include <boost/numeric/ublas/lu.hpp>
#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/progress.hpp>
#include <Eigen/Dense>
#include <iostream>
using namespace boost;
using namespace std;
const int n=9;
const int total=100000;
void test_ublas()
{
using namespace boost::numeric::ublas;
cout << "Boost.ublas ";
double r=1.0;
{
boost::progress_timer t;
for(int j=0;j!=total;++j)
{
//symmetric_matrix< double,lower,row_major,bounded_array<double,(1+n)*n/2> > A(n,n);
matrix<double,row_major,bounded_array<double,n*n> > A(n,n);
permutation_matrix< unsigned char,bounded_array<unsigned char,n> > P(n);
bounded_vector<double,n> v;
for(int i=0;i!=n;++i)
for(int k=0;k!=n;++k)
A(i,k)=0.0;
for(int i=0;i!=n;++i)
{
A(i,i)=1.0+i;
v[i]=i;
}
lu_factorize(A,P);
lu_substitute(A,P,v);
r+=inner_prod(v,v);
} …在Eigen中是否存在使用相对和绝对容差(numpy.allclose)比较矢量(矩阵)的函数？如果其中一个向量非常接近零,则标准isApprox失败.
在MATLAB中,您可以使用该函数计算矩阵的Jordan正规形式jordan.
它有NumPy和SciPy中的等效功能吗？
我正在计算Vandermonde matrix一个相当大的1D阵列.这样做的自然而干净的方法就是使用np.vander().但是,我发现这是约.比基于列表推导的方法慢2.5倍.
In [43]: x = np.arange(5000)
In [44]: N = 4
In [45]: %timeit np.vander(x, N, increasing=True)
155 µs ± 205 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
# one of the listed approaches from the documentation
In [46]: %timeit np.flip(np.column_stack([x**(N-1-i) for i in range(N)]), axis=1)
65.3 µs ± 235 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [47]: np.all(np.vander(x, N, …