q**j
发帖数: 10612 | 1 【 以下文字转载自 Statistics 讨论区 】
发信人: qqzj (小车车), 信区: Statistics
标 题: 求助：Matlab里面如何分解矩阵？
发信站: BBS 未名空间站 (Fri Jan 28 13:02:47 2011, 美东)
打算把一个positive definite, symmetric 的矩阵分解成两个矩阵， i.e.，
A = B*B.
按照原理可以用eigen decomposition 把A分成 V*D*V^(-1)。其中D是diagonal的。因为
A是正定的所以D的对角线上面都是正数。所以定义E是D的平方根，
A= V*D*V(-1) = V*E*E*V(-1) = V*E*V(-1) * V*E*V(-1)，所以定义
B = V*E*V(-1)，就可以有A = B*B。
可是在matlab里面用了 [V,D]=eig(A)以后，出来的D不是diagonal的。应该是计算精度
的问题吧。这下怎么办？有没有好办法来对付这个问题？多谢了。 |
t***s
发帖数: 4666 | 2 你是做什么的？为什么你的问题我都明白？hehe。eig by definition produce a
diagonal matrix. provide an exmaple of your matrix A?
因为
【在 q**j 的大作中提到】
: 【 以下文字转载自 Statistics 讨论区 】
: 发信人: qqzj (小车车), 信区: Statistics
: 标 题: 求助：Matlab里面如何分解矩阵？
: 发信站: BBS 未名空间站 (Fri Jan 28 13:02:47 2011, 美东)
: 打算把一个positive definite, symmetric 的矩阵分解成两个矩阵， i.e.，
: A = B*B.
: 按照原理可以用eigen decomposition 把A分成 V*D*V^(-1)。其中D是diagonal的。因为
: A是正定的所以D的对角线上面都是正数。所以定义E是D的平方根，
: A= V*D*V(-1) = V*E*E*V(-1) = V*E*V(-1) * V*E*V(-1)，所以定义
: B = V*E*V(-1)，就可以有A = B*B。
|
q**j
发帖数: 10612 | 3 正常情况下当然都是对角的了？关键是矩阵大了就不是这么理想了。我的矩阵都是3000
by 3000的。用了eig(A)以后给的对角矩阵不是完全对角的。而是block diagonal的。
每个block上面对角线上是正的实数加上一个很小的imaginary part，但是其他地方是
很小的虚数。我在想是否可以把
这些虚数硬变成0算了。
【在 t***s 的大作中提到】
: 你是做什么的？为什么你的问题我都明白？hehe。eig by definition produce a
: diagonal matrix. provide an exmaple of your matrix A?
:
: 因为
|
t***s
发帖数: 4666 | 4 try svd, it's more robust.
3000
【在 q**j 的大作中提到】
: 正常情况下当然都是对角的了？关键是矩阵大了就不是这么理想了。我的矩阵都是3000
: by 3000的。用了eig(A)以后给的对角矩阵不是完全对角的。而是block diagonal的。
: 每个block上面对角线上是正的实数加上一个很小的imaginary part，但是其他地方是
: 很小的虚数。我在想是否可以把
: 这些虚数硬变成0算了。
|
s*****v
发帖数: 360 | 5 Try SVD decomposition.
因为
【在 q**j 的大作中提到】
: 正常情况下当然都是对角的了？关键是矩阵大了就不是这么理想了。我的矩阵都是3000
: by 3000的。用了eig(A)以后给的对角矩阵不是完全对角的。而是block diagonal的。
: 每个block上面对角线上是正的实数加上一个很小的imaginary part，但是其他地方是
: 很小的虚数。我在想是否可以把
: 这些虚数硬变成0算了。
|
t***s
发帖数: 4666 | 6 or you can try sqrtm
【在 t***s 的大作中提到】
: try svd, it's more robust.
:
: 3000
|
o****e
发帖数: 4946 | 7 既然是数值计算,误差当然不可避免
你把D写成D ＝ E＊E + R， E是对角阵,R是误差
于是A ＝ V*E*E*V^(-1) + V*R*V^(-1)
令B=V*E*V^(-1),所以A=B*B + V*R*V^(-1)
最后一项是误差项, 算出来看看这样的误差是否可以接受就是了
3000
【在 q**j 的大作中提到】
: 正常情况下当然都是对角的了？关键是矩阵大了就不是这么理想了。我的矩阵都是3000
: by 3000的。用了eig(A)以后给的对角矩阵不是完全对角的。而是block diagonal的。
: 每个block上面对角线上是正的实数加上一个很小的imaginary part，但是其他地方是
: 很小的虚数。我在想是否可以把
: 这些虚数硬变成0算了。
|
q**j
发帖数: 10612 | 8 多谢，你这个办法好像和我的办法是一样的 :)
sqrtm很好很强大，虽然很费时间。估计matlab里面又在搞什么优化了。我还是不自己算
SVD了。出力不讨好的事情呀。
多谢三位指教。
【在 o****e 的大作中提到】
: 既然是数值计算,误差当然不可避免
: 你把D写成D ＝ E＊E + R， E是对角阵,R是误差
: 于是A ＝ V*E*E*V^(-1) + V*R*V^(-1)
: 令B=V*E*V^(-1),所以A=B*B + V*R*V^(-1)
: 最后一项是误差项, 算出来看看这样的误差是否可以接受就是了
:
: 3000
|
o****e
发帖数: 4946 | 9 sqrtm一样有误差,天下木有白吃的午餐
己算
【在 q**j 的大作中提到】
: 多谢，你这个办法好像和我的办法是一样的 :)
: sqrtm很好很强大，虽然很费时间。估计matlab里面又在搞什么优化了。我还是不自己算
: SVD了。出力不讨好的事情呀。
: 多谢三位指教。
|
q**j
发帖数: 10612 | 10 但是那个就是matlab的问题，不是我的了！
【在 o****e 的大作中提到】
: sqrtm一样有误差,天下木有白吃的午餐
:
: 己算
|
o****e
发帖数: 4946 | 11 ft，这个工作态度不professional
【在 q**j 的大作中提到】
: 但是那个就是matlab的问题，不是我的了！
|
M*******c
发帖数: 4371 | 12 Defined the smallest 实数 ep. Anything less than ep be ignored.
I remember there is a fuction called Re()
3000
【在 q**j 的大作中提到】
: 正常情况下当然都是对角的了？关键是矩阵大了就不是这么理想了。我的矩阵都是3000
: by 3000的。用了eig(A)以后给的对角矩阵不是完全对角的。而是block diagonal的。
: 每个block上面对角线上是正的实数加上一个很小的imaginary part，但是其他地方是
: 很小的虚数。我在想是否可以把
: 这些虚数硬变成0算了。
|
