wwdjkp 发表于 2006-8-22 16:47

关于 matlab 里的矩阵除法问题

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

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

xinyuxf 发表于 2006-8-23 10:13

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

happy 发表于 2006-8-23 15:22

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

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

wwdjkp 发表于 2006-8-23 16:43

原帖由 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, 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 发表于 2006-8-23 18:06

I'm not sure...It seems like this???
b=round(10*rand(3));
A=round(10*rand(3,4));
X=pinv(A)*b

happy 发表于 2006-8-23 19:10

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

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

举一个简单的例子,比如



a=;b=;


则采用两种方法求解,x=a\b 和 x=pinv(a)b

第一种方法的结果是:
第二种方法的结果是:

happy 发表于 2006-8-23 19:13

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

wwdjkp 发表于 2006-8-24 18:34

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

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


的确如此,比如取M=325,求解出来的X大多数行都是0,如以下截取的前13行,为表示方便,转置了一下
0 0-1.2593e+0090 0 0 0 0 0 0 0 0 -1.3489e+009
0 0-6.8914e+0090 0 0 0 0 0 0 0 0 -7.1698e+009
0 0-1.7833e+0100 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
=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 发表于 2006-8-24 20:53

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

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

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

thierryhenry 发表于 2008-3-25 18:59

受益匪浅

受益匪浅:handshake
页: [1]
查看完整版本: 关于 matlab 里的矩阵除法问题