[求助]关于用YACOBI旋转方法来求解下面问题

allen

n=7;
m_trz=[1 3; 1 4; 4 6]; % 元素不为零的下标位置；

%P为参考矩阵，就是说M最后会形成这样的形式；
%意思是下面的M(1,3),M(1,4)，M(4,6)不为零。消去M中其它对应P矩阵为0的元素。
P＝[1 1 1 0 0 0 0 0 0;
       1 1 1 0 0 0 0 0 0;
       1 1 1 1 1 0 0 0 0;
       0 0 1 1 1 0 0 0 0;
       0 0 1 1 1 1 0 0 0;
       0 0 0 0 1 1 1 0 0;
       0 0 0 0 0 1 1 1 0;
       0 0 0 0 0 0 1 1 1;
       0 0 0 0 0 0 0 1 1];


M_initial=[ 0.0211 0.8886 -0.0000  -0.0000  0.0000  0.0000     0;
                  0.8886 0.0252  0.6162   0.0000 -0.0000  0.0942     0;
                -0.0000  0.6162  0.0187  0.5103   0.1880  0.0700   0.0000;
                -0.0000  0.0000  0.5103 -0.4857  0.4546            0   0.0000;
                 0.0000 -0.0000  0.1880  0.4546 -0.0228  0.6115  -0.0000;
                 0.0000  0.0942   0.0700          0   0.6115  0.0267   0.8881;
                          0           0   0.0000 0.0000  -0.0000  0.8881   0.0211];

利用YACOBI旋转方法，得到和P矩阵类似的矩阵，在YACOBI方法中，要输入一串的theta值，这为优化变量，
请问怎么利用MATLAB的优化算法来实现这个要求？
要求进行YACOBI变换前后的M矩阵的特性不变。M为关于主线对称的方阵。

happy

没看懂你在说什么

allen

晕。就是想建立一个优化算法，输入不同的THETA，等到所有P为0的元素对应M矩阵的值的和的平方为零时（或者小于一个10E-8，这个是建立的评价函数），。输出这个最后的M矩阵。

就是这个意思。可能相关的，要看JAKOBI旋转法。请查看附件。

但里面只是要得到一个对称的只有主对角线存有元素，和我的期望得到的结果不太一样。

happy

看到你的短信了，你到底是要求最小呢？还是说u_cost<10e-8就行？

happy

但里面只是要得到一个对称的只有主对角线存有元素

这里的里面指的是什么？看来我的理解能力确实太差了点

allen

不好意思。我要求的是u_cost<10e-8.

里面是指附件里面的例子计算的结果。

但是我需要的结果是，M矩阵类似如P0矩阵一样，当P0＝1，M为大于10e-8；当P0＝0，M为小于10e-8。

不知道说清楚没，谢谢版主！

allen

版主，能帮我解答一下么，谢谢了！

或者，哪位线性代数高手帮我看一下啊，谢过大家了。拜托

bainhome

JACOBI法严格意义上讲是线性代数理论在数值分析中特征值特征向量计算里的一个应用，不完全是线性代数的问题。
你的附件里提到的是基本理论，实际计算根本不需要用什么优化算法去求该θ值，有其自己的计算方式，所以我说你看书只看了一半，或者没有自己查阅资料了解问题全貌就匆忙发问！
由于无论是多少阶的矩阵，其每次迭代的变换旋转矩阵中只有两行不同，所以是一个关于tanθ的二次方程，由jacobi法的收敛性可以证明进行一系列旋转变换最终理论上趋近于一个对角阵，该对角阵的主对角线上元素即为特征值，根据该计算方法，我简单写如下程序：

1. function JacobiSolveLambda(A,t)
   
2. % A必须为一个对称矩阵
   
3. % t:迭代终止控制误差
   
4. clc
   
5. count=0;
   
6. A1=triu(A,1);
   
7. T1=triu(A,1)+tril(A,-1);
   
8. T=norm(T1(:));
   
9. while T>t
   
10.     A2=abs(A1);
    
11.     A2=A2(:);
    
12.     maxA2=max(A2);
    
13.     A2=reshape(A2,size(A,1),size(A,2));
    
14.     [i1,j1]=find(A2==maxA2);
    
15.     if length(i1)>1
    
16.         i1=i1(1);j1=j1(1);
    
17.     end
    
18.     d=(A(i1,i1)-A(j1,j1))/(2*A(i1,j1));
    
19.     if d>=0
    
20.         sd=1;
    
21.     else
    
22.         sd=-1;
    
23.     end
    
24.     tt=sd/(abs(d)+sqrt(d^2+1));
    
25.     c=1/sqrt(1+tt^2);
    
26.     s=c*tt;
    
27.     P=eye(length(A));
    
28.     P(i1,i1)=c;
    
29.     P(j1,j1)=c;
    
30.     P(i1,j1)=s;
    
31.     P(j1,i1)=-s;
    
32.     P;
    
33.     A=P*A*P';
    
34.     A1=triu(A,1);
    
35.     count=count+1;
    
36.     T1=triu(A,1)+tril(A,-1);
    
37.     T=norm(T1(:));
    
38. end
    
39. disp('特征值:')
    
40. S=diag(A)
    
41. disp(['迭代次数n=',num2str(count)])

复制代码

command windows中运行：

1. A=round(10*rand(4))
   
2. A1=triu(A,1)
   
3. A=diag(diag(A))+A1+A1'
   
4. JacobiSolveLambda(A,1e-2)
   
5. eig(A)

复制代码

可以看出和eig命令得到的结果基本相同。

allen

谢谢bainhome的回复！很是受用。
但，我要求的结果不是得到特征值，我要求的矩阵是
当P为1时，M矩阵在应用旋转法后保留有大于控制误差的值；
当P为0时，M矩阵在应用旋转法后不保留大于控制误差的值，无限趋近于0。

这是我要求的结果，而且，M矩阵关于对角线对称的，得到的结果也必须是关于对角线对称的。
也就是M(1,2)=M(2,1).....

不好意思，麻烦您了。因为我实在是对代数这方面不是很精通，但是工程应用需要，所以只要发帖求助了。

bainhome

我也不是学数学的，不过又反复看了几遍你的问题，感觉不能帮你。
原先以为你想编写特征值求解的程序，只是不知道旋转变换如何实现，想通过求极值的方法得到aij=aji=0时的θ值。而原贴又没有说明和强调这一点。以至于有上面我的帖子的误解。恕我直言，第一个帖子根本没问得突出到重点上：
对矩阵做JACOBI旋转变换只是把对称阵转化为对角阵而已，你附件里的内容基本就是数值分析中求特征值的JACOBI法的原话，因此在进行旋转变换的一系列过程中M矩阵肯定是对称的，用我的程序随便试，每次迭代中都是对称阵
这是你的帖子原文：

n=7;
m_trz=[1 3; 1 4; 4 6]; % 元素不为零的下标位置；

%P为参考矩阵，就是说M最后会形成这样的形式；
%意思是下面的M(1,3),M(1,4)，M(4,6)不为零。消去M中其它对应P矩阵为0的元素。
P＝[9×9];
M_initial=[ 7×7];
利用YACOBI旋转方法，得到和P矩阵类似的矩阵，在YACOBI方法中，要输入一串的theta值，这为优化变量，
请问怎么利用MATLAB的优化算法来实现这个要求？
要求进行YACOBI变换前后的M矩阵的特性不变。M为关于主线对称的方阵。

1.P和M一个7阶一个9阶，维数不同如何旋转？
2."意思是下面的M(1,3),M(1,4)，M(4,6)不为零。消去M中其它对应P矩阵为0的元素"，里面“消去”是什么意思？缩维还是置零？“M(1,3),M(1,4)，M(4,6)不为零”，我看M里非对角元素不为零的不止这几项，不知什么意思？
3.“在YACOBI方法中，要输入一串的theta值”。
这里的“一串θ值”未必准确吧？？？这个问题的求法就在我的程序里，每次迭代时只出现一个θ值。而且是直接求cos(θ)和sin(θ)的数值，JACOBI旋转法在计算机实现的时候我理解的是根本不需要去直接求这个θ值，它的一系列迭代是通过“P1P2P3...PnAPn'...P3'P2'P1'”来实现的，角度的正、余弦算出之后，P就自然得到，剩下就是矩阵乘了，因为JACOBI的基本思想就是通过一组平面旋转变换也就是正交相似变换将对称阵转化为对角阵，恕我愚钝，从我理解，没看出来jacobi旋转变换中有对θ角的极值求解和优化的意思。充其量是个以tanθ为变量的方程，而且是算法中直接套用的公式。并未用solve去求解，因为每回迭代都是那个方程，提前解完了事儿！
4.“当P为1时，M矩阵在应用旋转法后保留有大于控制误差的值；
当P为0时，M矩阵在应用旋转法后不保留大于控制误差的值，无限趋近于0。”
这是你最后的帖子里的内容，我越看越糊涂！
4-1.P是个矩阵，什么叫P为1或者为0，行列式还是每个元素？为什么不写清楚！
4-2.“保留有大于控制误差的值”，如何保留？我猜测是认为矩阵M旋转变换后得到的Mn矩阵里，大于迭代精度的元素不变，小于迭代精度的元素置0，如果是这意思，为什么要写得不明不白，如果不是这个意思，请教下你到底要表达什么？
4-3.“无限趋近于0”是什么意思，在数值计算过程中我没见过所谓“趋近于0”的表达方式，只有在某个阈值以下认为是0，以上不为0
5.最初的帖子里你说经过变换后的M要形成一个类似P的矩阵形式，而且是对称阵，你自己看看，你给的样板阵P，其本身根本就不是一个对称阵，更何况和所给M不同维！
我的建议是你自己先消化一下你的内容，或者举一个简单的例子把你认为要优化的，关于θ的表达式萃取出来再与其他人讨论，如果是一个优化问题，可以看看fmincon、fminbnd、fminunc三个函数，如果表达式固定且容易求导，可通过diff命令得到其导函数的表达形式，解一个关于θ的非线性方程也能得到这个值。
不过不客气地说，你让人看这么多资料，在网上能得到一个可用结果的概率基本是零。而且除非有人做过完全相同的问题，否则以你表达的内容来说，大多数还不知道你想干什么，既然不是特征值的问题，那我也是一头雾水，实在无能为力，抱歉。

[ 本帖最后由 bainhome 于 2006-8-20 04:09 编辑 ]

allen

老大，呵呵，感谢你的回问，在这里受教了。我自己有点急，呵呵，所以问的问题错误一堆，没去审核，实在是报歉。

1，不好意思，这个是我写错！P和M应该是同维矩阵；P只是一个参考，这里的P应该是7*7维的。
2，我不知道JACOBI旋转方法是否可以达到置零，应该是小于迭代精度的元素吧。
　最后想得到的M矩阵，应该是一个三对角矩阵，另外再加上这三个元素M(1,3),M(1,4)，M(4,6)不为零，而其它元素都应小于迭代精度元素。

3,其实我现在想要的，只是利用JACOBI旋转的特性，对矩阵M进行处理，但是处理前后它的特征值及矩阵特性不变，这里看过资料，说要对θ值进行优化来得到。所以直接提出来；
4-1 是元素；
4-2 是想通过旋转方法来消去问题2中不想要的位置的元素。
4-3 应该是在某个阈值以下认为是0，以上不为0

5  P矩阵应该是：
P＝ 1     1     1     1     0     0     0
     1     1     1     0     0     0     0
     1     1     1     1     0     0     0
     1     0     1     1     1     1     0
     0     0     0     1     1     1     0
     0     0     0     1     1     1     1
     0     0     0     0     0     1     1

真是感谢你的这些建议，可能我是理解比较含糊就帖出来了，我的问题是太多了，因为时间也太急，还请谅解。
请问有邮箱吗，我能和你联系吗？谢谢！欢迎到我论坛做客。呵呵，在这里有做广告之嫌，
还请版主不要责备。