关于 matlab 里的矩阵除法问题

wwdjkp

这个问题简单来说就是求解不定或超定方程A*X＝B，A为m*n的矩阵
在matlab里直接用X＝B\A就可以了
可我现在需要用C语言实现
所以请问大家matlab中是如何实现B\A的
或者在matlab 里的“ / ”除运算符调用的是哪个内部函数？
小弟不胜感激:@)

[ 本帖最后由 eight 于 2008-3-25 19:50 编辑 ]

xinyuxf

解这样的方程，是不是要用到数值计算中的迭代方法阿？

happy

在matlab中矩阵的左除"\"是一种自适应算法
对于有唯一解的线性方程组求出精确解
对于超定方程线性方程组求出最小二乘解
对于欠定线性方程组求其基本解

这个好像是个内部函数，看不了它的代码

wwdjkp

原帖由 happy 于 2006-8-23 15:22 发表
在matlab中矩阵的左除"\"是一种自适应算法
对于有唯一解的线性方程组求出精确解
对于超定方程线性方程组求出最小二乘解
对于欠定线性方程组求其基本解

这个好像是个内部函数，看不了它的代码

首选感谢教授的热心指点:@)

我现在面临的问题是欠定线性方程组的求解，即A为M*(M+1)的长方阵，X为(M+1)*3的矩阵，B为M*3的矩阵
请问教授在Matlab 里是如何求此基本解的？

我最近也看了不少广义逆矩阵的资料，有A-  A+  A(1,3)  A(1,4)等等，拭了几种都跟matlab 算的不一样，另外还有用QR分解的方法，现在正在做，不知道行不行

下面是matlab 里的mldivide文件，不知道是否是调用的这个函数来求"\"运行符，请指教

function X = mldivide(A, B)
%MLDIVIDE Symbolic matrix left division.
%   MLDIVIDE(A,B) overloads symbolic A \ B.
%   X = A\B solves the symbolic linear equations A*X = B.
%   Warning messages are produced if X does not exist or is not unique.  
%   Rectangular matrices A are allowed, but the equations must be
%   consistent; a least squares solution is not computed.
   
%   Copyright 1993-2003 The MathWorks, Inc.
%   $Revision: 1.18.4.2 $  $Date: 2004/04/16 22:22:54 $

A = sym(A);
B = sym(B);

if all(size(A) == 1)
   % Division by a scalar
   X = ldivide(A,B);

elseif ndims(A) > 2 | ndims(B) > 2
   error('symbolic:sym:mldivide:errmsg1','Input arguments must be 2-D.')

elseif size(A,1) ~= size(B,1)
   error('symbolic:sym:mldivide:errmsg2','First dimensions must agree.')

else
   % Matrix divided by matrix

   X = maple('linsolve',char(A),char(B),'''_rank''');

   % Solution does not exist.
   if isempty(X)
      warning('symbolic:sym:mldivide:warnmsg1','System is inconsistent. Solution does not exist.')
      X = Inf;
      X = sym(X(ones(size(A,2),size(B,2))));
      maple('_rank := ''_rank'';');
      return
   end;

   % Check rank and clear _rank in Maple workspace.
   if str2double(maple('_rank')) < min(size(A))
      warning('symbolic:sym:mldivide:warnmsg2','System is rank deficient. Solution is not unique.')
   end
   maple('_rank := ''_rank'';');

   % Set any free parameters, _t[k][j], to zero.
   t = findstr(X,'_t[');
   s = findstr(X,']');
   for k = fliplr(t)
      r = s(s > k);
      X(k:r(2)) = '0';
   end
   X = maple('',sym(X));
end

bainhome

I'm not sure...It seems like this???

1. b=round(10*rand(3));
   
2. A=round(10*rand(3,4));
   
3. X=pinv(A)*b

复制代码

happy

对于欠定方程组matlab一般采用两种方法求解

1. 采用除法求解，这种方法的求解结果是具有最多零元素的解
2. 基于伪逆pinv的求解方法，这种方法求解的具有最小长度或范数的解，楼上所说的属于这种情况

举一个简单的例子，比如

a=[1 2 3;2 3 4];b=[1;2];

则采用两种方法求解，x=a\b 和 x=pinv(a)&#61482;b

第一种方法的结果是：[1，0，0]
第二种方法的结果是：[0.83，0.33，-0.17]

happy

mldivide这个函数一直没有注意过，不过从说明上看好像这个函数就是左除

wwdjkp

原帖由 happy 于 2006-8-23 19:10 发表
对于欠定方程组matlab一般采用两种方法求解

1. 采用除法求解，这种方法的求解结果是具有最多零元素的解

的确如此，比如取M=325，求解出来的X大多数行都是0，如以下截取的前13行，为表示方便，转置了一下
0 0  -1.2593e+009  0 0 0 0 0 0 0 0 0 -1.3489e+009
0 0  -6.8914e+009  0 0 0 0 0 0 0 0 0 -7.1698e+009
0 0  -1.7833e+010  0 0 0 0 0 0 0 0 0 -1.7629e+010

请教授大概讲一下计算原理，我得用C语言来实现

另外，今天我拭验了一下，似乎是采取如下的方法，下面是一段例程，可以直接运行，拷贝到matlab里看着清楚一些

clc
clear
m=15;
a=rand(m,m+1);
b=rand(m,3);
%相当于求解欠定线性方程组a*xx=b
[q r]=qr(a);%对a进行QR分解,即a=q*r,q为m*m的正交阵(重要性质：q*q'=I)，r为m*(m+1)的上三角矩阵
rr=r;%保存r
for k=1:m+1
    r=rr;
    r(:,k)=[];%去掉r阵的第k列，r成为方阵，从而可以用inv求逆
    X{1,k}=inv(r)*q'*b;
end
xx=a\b;%此为由matlab内部函数求出的方程解
%说明
%若xx矩阵中k行全为0，即表明去掉了r阵的k行，亦即 X{1,k}中的矩阵即为xx　请教授明查
%
%所以现在的问题就是不知道matlab是依据什么规则去掉r阵的某一行，

再次感谢教授的指导:@)
俺是研二学生，现在课题马上就要作完了，就卡在这一个问题上了，恳请教授不吝赐教

[ 本帖最后由 wwdjkp 于 2006-8-24 18:55 编辑 ]

happy

今天又查了一下相关资料，可以确定的是mldivide就是左除，右除是mrdivide

其基本算法采用的QR因式分解，超定方程的整个求解方法你可以查看qr这个命令的帮助文件Example 2.

欠定的类似的，你看一下相关数值分析的书籍吧

thierryhenry

受益匪浅:handshake