julia 令人发指的 Cholesky 分解

正在Debug，惊是不断，喜就免了。先来看看华丽的接口吧，一般人还找不到

x=rand(3,1);
m=x*x';
m=m+.001eye(3);
u=cholfact(m,:U)[:U].data;

这样才算是取到你想要的输出，功能够丰富吧；
这样的优雅，结果必定亮瞎狗眼：

julia> x=rand(3,1);m=x*x';m=m+.001eye(3);u=cholfact(m,:U)[:U].data;u'*u-m
3×3 Array{Float64,2}:
 0.014269    0.0208616   0.00302363
 0.0208616   0.00687607  0.00264472
 0.00302363  0.00264472  0.0       
julia> x=rand(3,1);m=x*x';m=m+.001eye(3);u=cholfact(m,:U)[:U].data;u'*u-m
3×3 Array{Float64,2}:
 0.245597    0.00818602   0.0347094 
 0.00818602  0.000256138  0.00112054
 0.0347094   0.00112054   0.0       

果然哭瞎在茅房。。。orz，　这次是真跪了，MD，程序都写完了，你给我看这个。
看看人家octave也比你这渣强到不知哪去了：

octave:7> x=rand(3,1);m=x*x';m=m+.001*eye(3);l=chol(m);l'*l-m
ans =

  -3.4694e-18   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00
octave:8> x=rand(3,1);m=x*x';m=m+.001*eye(3);l=chol(m);l'*l-m
ans =

   0.0000e+00  -1.1102e-16   0.0000e+00
  -1.1102e-16   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   0.0000e+00
octave:9> x=rand(3,1);m=x*x';m=m+.001*eye(3);l=chol(m);l'*l-m
ans =

   2.1684e-19   0.0000e+00   0.0000e+00
   0.0000e+00   0.0000e+00   1.7347e-18
   0.0000e+00   1.7347e-18   2.1684e-19

已经摊上这货了，怎么办？想到了numpy，只有这样试试了:

julia> using PyCall
julia> @pyimport numpy as np 
julia> npchol=np.linalg[:cholesky]     # 瞧瞧这华丽的修饰
julia> x=rand(3,1);m=x*x';m=m+.001eye(3);l=npchol(m);np.dot(l,l')-m
3×3 Array{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0

julia> x=rand(3,1);m=x*x';m=m+.001eye(3);l=npchol(m);np.dot(l,l')-m
3×3 Array{Float64,2}:
 -2.77556e-17  0.0  -3.46945e-18
  0.0          0.0   0.0        
 -3.46945e-18  0.0   0.0     

# 刚才被感动昏了，这样输入也是可以的
julia> x=rand(3,1);m=x*x';m=m+.001eye(3);l=npchol(m);l*l'-m
3×3 Array{Float64,2}:
 2.1684e-19  0.0          0.0
 0.0         5.55112e-17  0.0
 0.0         0.0          0.0

做完这些，终于依稀想起上期做LU　factorization的动人情形了，不想验证了，心累。

转载于:https://www.cnblogs.com/chenyliang/p/6817497.html