标签: linear-algebra

我想反转一个4x4矩阵.我的数字以定点格式存储(准确地说是1.15.16).

使用浮点运算,我通常只需构建伴随矩阵并除以行列式(例如强力解决方案).到目前为止,这对我有用,但在处理定点数时,由于使用了所有乘法,我得到了不可接受的精度损失.

注意:在定点运算中,我总是丢弃一些立即结果的最低有效位.

那么 - 反转矩阵的数值最稳定的方法是什么？我不太关心性能,但简单地说浮点就是放慢我的目标架构.

我正在用Python原型化一个新系统; 功能主要是数字.

一个重要的要求是能够使用不同的线性代数后端:从单个用户实现到通用库,如Numpy.线性代数实现(即后端)必须独立于接口.

我最初的架构尝试如下:

(1)定义系统界面

>>> v1 = Vector([1,2,3])
>>> v2 = Vector([4,5,6])
>>> print v1 * v2
>>> # prints "Vector([4, 10, 18])"

(2)实现允许独立于后端使用该接口的代码

# this example uses numpy as the back-end, but I mean
# to do this for a general back-end
import numpy 
def numpy_array(*args): # creates a numpy array from the arguments
    return numpy.array(*args)

class VectorBase(type):
    def __init__(cls, name, bases, attrs):
        engine = attrs.pop("engine", None)
        if not engine:
            raise RuntimeError("you need to specify an engine")
        # …

我试图将包含\运算符的代码从Matlab(Octave)转换为Python.示例代码

B = [2;4]
b = [4;4]
B \ b

这有效并产生1.2作为答案.使用此网页

http://mathesaurus.sourceforge.net/matlab-numpy.html

我翻译为:

import numpy as np
import numpy.linalg as lin
B = np.array([[2],[4]])
b = np.array([[4],[4]])
print lin.solve(B,b)

这给了我一个错误:

numpy.linalg.linalg.LinAlgError: Array must be square

为什么Matlab \与B的非方矩阵一起工作？

对此有何解决方案？

我在Matlab上使用了一个函数:

[V,D] = eig(C);

我看到V和D总是按升序排序.它是否总是那样,或者我应该在获得V和D价值后对它们进行排序？

我有两个图像,并使用筛选找到三个相似的2D点.我需要计算图像之间的仿射变换.不幸的是,我错过了讲座,那里的信息对我来说有点密集.计算这个2x3矩阵的一般方法是什么？

我有一个2x3矩阵中的点矩阵[x1 y1; x2 y2; x3 y3],但我从那里迷失了.谢谢你的帮助.

在MATLAB中,我使用归一化的八点算法计算了基本矩阵(两个图像).由此我需要对3D空间中的相应图像点进行三角测量.根据我的理解,要做到这一点,我需要旋转和翻译图像的相机.最简单的方法当然是首先校准相机然后拍摄图像,但这对我的应用来说太紧张了,因为它需要这个额外的步骤.

这样我就可以进行自动(自拍)相机校准.我看到了捆绑调整的提及,但是在"邀请3D视觉"中它似乎需要初始平移和旋转,这让我觉得需要校准相机或者我的理解不足.

所以我的问题是如何自动提取旋转/平移,以便我可以将图像点重新投影/三角测量到3D空间.任何MATLAB代码或伪代码都会很棒.

我是大学的CS专业,为我的Calc III课程编写一个涉及奇异值分解的编程项目.该想法基本上是将m×n维度的图像转换为m×n矩阵,其中每个元素是表示点(m,n)处的像素的颜色通道(r,g,b)的元组.我正在使用Python,因为它是我迄今为止(好)教授的唯一语言.

据我所知,Python通常不喜欢元组作为数组的元素.我对自己进行了一些研究并找到了一种解决方法,即按如下方式预先分配数组:

def image_to_array(): #converts an image to an array  
    aPic = loadPicture("zorak_color.gif")  
    ph = getHeight(aPic)  
    pw = getWidth(aPic)  
    anArray = zeros((ph,pw), dtype='O')  
    for h in range(ph):  
         for w in range(pw):             
            p = getPixel(aPic, w, h)  
            anArray[h][w] = (getRGB(p))  
    return anArray

这适用于赋值的第一部分,它只是将图像转换为矩阵(不涉及线性代数).

然而,具有SVD的部分变得更加棘手.当我调用内置的numPy svd函数时,使用我从我的图像构建的数组(其中每个元素是一个元组),我得到以下错误:

Traceback (most recent call last):
  File "<pyshell#5>", line 1, in -toplevel-
    svd(x)
  File "C:\Python24\Lib\site-packages\numpy\linalg\linalg.py", line 724, in svd
    a = _fastCopyAndTranspose(t, a)
  File "C:\Python24\Lib\site-packages\numpy\linalg\linalg.py", line 107, in _fastCopyAndTranspose
    cast_arrays = cast_arrays + (_fastCT(a.astype(type)),)
ValueError: setting …

由于D很大,我没有足够的内存来简单地创建对角D-by-D矩阵.我不断收到"内存不足"的错误.

我没有在第一次乘法中执行M x D x D操作,而是执行M x D操作,但我的代码仍然需要很长时间才能运行.

任何人都可以找到一种更有效的方法来执行乘法运算A'*B*A吗？这是我到目前为止所做的尝试:

D=20000
M=25

A = floor(rand(D,M)*10);
B = floor(rand(1,D)*10);

for i=1:D
    for j=1:M
        result(i,j) = A(i,j) * B(1,j);
    end
end    

manual = result * A';
auto = A*diag(B)*A';
isequal(manual,auto)

我在解决Ax = B形式的系统时遇到了麻烦

该系统的解决方案应该是

x = inv(A)*B

但是,这不起作用.

我尝试上面的代码行时收到以下错误消息:

Warning: Matrix is close to singular or badly scaled.
     Results may be inaccurate. RCOND = 1.156482e-018. 

似乎matlab在反转我指定的矩阵时遇到了麻烦.我试着通过输入inv(A)*A来验证反函数是否正常工作

这应该给出单位矩阵,但是我得到了相同的错误和一些垃圾数.

这是我正在使用的A矩阵:

A = [5/2   1/2  -1     0     0    -1/2  -1/2   0     0
     1/2   1/2   0     0     0    -1/2  -1/2   0     0 
    -1     0     5/2  -1/2  -1     0     0    -1/2   1/2
     0     0    -1/2   1/2   0     0     0     1/2  -1/2
     0     0    -1     0     3/2  -1/2   1/2   0     0
    -1/2  -1/2   0     0    -1/2 …

我在R中有一个计算大量最小二乘解决方案(> 10,000:通常为100,000+)的程序,并且在分析之后,这些是该程序的当前瓶颈.我有一个矩阵A,其列向量对应于生成向量和解决方案b.我试图求解最小二乘解x的Ax=b.矩阵的大小通常为4xj - 其中很多都不是正方形(j <4),因此对于欠定系统的一般解决方案正是我所寻求的.

主要问题:在R中解决欠定系统的最快方法是什么？我有很多利用Normal Equation的解决方案,但我在R中寻找比下面任何一种方法更快的例程.

例如:x通过Ax = b给定以下约束来解决给定的系统:

• 系统没有必要确定[通常不确定](ncol(A)<= length(b)始终保持).因此ncol (A) <= length(b)无效,因为求解需要方阵.
• 您可以假设solve(A,b)(等效于t(A) %*% A)是非单数的 - 在程序的早期检查它
• 您可以使用R中免费提供的任何包装
• 解决方案不一定非常漂亮 - 只需要快速
• 大小上限crossprod(A)合理地为10x10,零元素很少发生 - A通常非常密集

两个随机矩阵用于测试......

A = matrix(runif(12), nrow = 4)
b = matrix(runif(4), nrow = 4)

以下所有功能均已分析.它们转载于此:

f1 = function(A,b)
{
  solve(t(A) %*% A, t(A) %*% b)
}
f2 = function(A,b)
{
  solve(crossprod(A), …