{"id":5300,"date":"2014-01-09T09:08:41","date_gmt":"2014-01-09T00:08:41","guid":{"rendered":"http:\/\/www.moonmile.net\/blog\/?p=5300"},"modified":"2014-01-07T14:10:15","modified_gmt":"2014-01-07T05:10:15","slug":"f-lu%e5%88%86%e8%a7%a3%e3%82%92%e4%bd%9c%e3%82%8b","status":"publish","type":"post","link":"http:\/\/www.moonmile.net\/blog\/archives\/5300","title":{"rendered":"[F#] LU\u5206\u89e3\u3092\u4f5c\u308b"},"content":{"rendered":"

\n\u65e7\u30bf\u30a4\u30d7\u306a\u306e\u3067\u3001\u9023\u7acb\u4e00\u6b21\u65b9\u7a0b\u5f0f\u3092\u89e3\u304f\u3068\u304d\u306b\u9006\u884c\u5217\u3092\u4f7f\u3063\u3066\u89e3\u3051\u3070OK…\u3068\u982d\u304b\u3089\u601d\u3044\u8fbc\u3093\u3067\u3044\u305f\u306e\u3067\u3059\u304c\u3001Gauss\u306e\u6d88\u53bb\u6cd5\u3092\u4f7f\u3046\u3068\u9006\u884c\u5217\u306f\u51fa\u3066\u3053\u306a\u304f\u3066\u300c\uff1f\u300d\u3068\u306a\u308b\u308f\u3051\u3067\u3059\u3002\u4e0a\u4e09\u89d2\u5316\u3092\u4f7f\u3046\u306e\u3060\u304b\u3089\u3001LU\u5206\u89e3\u3067\u3044\u3044\u3093\u3067\u3059\u3088\u306d\u3002\n<\/p>\n

\nLU\u5206\u89e3\u3067\u884c\u5217\u5f0f\u3068\u9006\u884c\u5217\u306e\u8a08\u7b97\u3000\u305d\u306e\uff13: \u30e1\u30e2\u30d6\u30ed\u30b0
\nhttp:\/\/sssiii.seesaa.net\/article\/379560261.html<\/a>\n<\/p>\n

\n\u306e\u3068\u3053\u308d\u304b\u3089\u3001\n<\/p>\n

\n\u9023\u7acb\uff11\u6b21\u65b9\u7a0b\u5f0f I
\nhttp:\/\/nalab.mind.meiji.ac.jp\/~mk\/labo\/text\/linear-eq-1.pdf<\/a>\n<\/p>\n

\n\u306ePDF\u3092\u8aad\u3093\u3067\u3001\u306a\u308b\u307b\u3069\u3001LU\u5206\u89e3\u3067\u5341\u5206\u3067\u3059\u306d\u3002\u3063\u3066\u3053\u3068\u3067\u81ea\u524d\u3067\u5b9f\u88c5\u3057\u3066\u307f\u307e\u3059\u3002\n<\/p>\n

\n4 LU\u5206\u89e3
\nhttp:\/\/akita-nct.jp\/yamamoto\/lecture\/2004\/5E\/linear_equations\/text\/html\/node4.html#eq:LU____________<\/a>\n<\/p>\n

\n\u3092\u53c2\u8003\u306b\u3057\u3066\u5f0f\u3092\u6d41\u7528\u3057\u307e\u3059\u3002\u6700\u521d\u3001\u624b\u9806\u304c\u3046\u307e\u304f\u308f\u304b\u3089\u306a\u304f\u3066\u3001\u81ea\u5206\u3067 4×4 \u306e\u884c\u5217\u3092\u4f5c\u3063\u3066\u624b\u8a08\u7b97\u3092\u3057\u305f\u5f8c\u306b\u3001\u6570\u5f0f\u306b\u76f4\u3057\u3066\u3001\u5143\u306e\u5f0f\u3068\u4e00\u81f4\u3059\u308b\u3053\u3068\u3092\u78ba\u8a8d\u3002\u306a\u308b\u307b\u3069\u3001\u884c\u5358\u4f4d\u3058\u3083\u306a\u304f\u3066\u300c\u5217\u5358\u4f4d\u300d\u3067\u8a08\u7b97\u3059\u308b\u306e\u3067\u3001\u8a08\u7b97\u3059\u308b\u9806\u756a\u306f\u3001\u03b211,\u03b121,\u03b131,\u03b141 \u3067\u8a08\u7b97\u3057\u305f\u5f8c\u306b\u3001\u03b212,\u03b222,\u03b132,\u03b142 \u3068\u8a08\u7b97\u3057\u3066\u3044\u304f\u3093\u3067\u3059\u306d\u3002\u3082\u3046\u4e00\u5ea6\u8aad\u307f\u76f4\u3057\u305f\u3089\u305d\u3046\u66f8\u3044\u3066\u3042\u308b\u3057…\n<\/p>\n

\n\u51fa\u6765\u4e0a\u304c\u3063\u305f\u30b3\u30fc\u30c9\u304c\u3053\u3093\u306a\u611f\u3058\u3002\n<\/p>\n

[code lang=”csharp”]
\n\/\/\/ LU\u5206\u89e3
\nlet LU ( A : matrix ) =
\n let n = A.NumCols
\n let L = Matrix.identity n
\n let U = Matrix.zero n n<\/p>\n

for j in 0..n-1 do
\n \/\/ U\u3092\u8a08\u7b97
\n for i in 0..j do
\n U.[i,j] <- A.[i,j]\n for k in 0..i-1 do\n U.[i,j] <- U.[i,j] - L.[i,k]*U.[k,j]\n \/\/ L\u3092\u8a08\u7b97\n for i in j+1..n-1 do\n L.[i,j] <- A.[i,j]\n for k in 0..j-1 do\n L.[i,j] <- L.[i,j] - L.[i,k]*U.[k,j]\n L.[i,j] <- L.[i,j] \/ U.[j,j]\n \/\/ \u7d50\u679c\u3092\u8fd4\u3059\n (L,U)\n[\/code]\n\n\n

\n\u6b21\u306e\u3088\u3046\u306b\u884c\u5217\u3092\u8a2d\u5b9a\u3057\u3066\u304a\u304f\u3068\n<\/p>\n

[code lang=”csharp”]
\nlet A = matrix [[1.0;2.0;1.0];
\n [2.0;1.0;0.0];
\n [1.0;1.0;2.0]]<\/p>\n

LU A
\n[\/code]<\/p>\n

\n\u3053\u3093\u306a\u611f\u3058\u3067\u3001L \u3068 U \u306b\u3057\u3066\u8fd4\u3057\u3066\u304f\u308c\u307e\u3059\u3002\n<\/p>\n

[code lang=”csharp”]
\nval it : matrix * matrix =
\n (matrix [[1.0; 0.0; 0.0]
\n [2.0; 1.0; 0.0]
\n [1.0; 0.3333333333; 1.0]], matrix [[1.0; 2.0; 1.0]
\n [0.0; -3.0; -2.0]
\n [0.0; 0.0; 1.666666667]])
\n[\/code]<\/p>\n

\n\u5148\u306e\u30da\u30fc\u30b8\u306b\u3082\u66f8\u3044\u3066\u3042\u308b\u306e\u3067\u3059\u304c\u3001Ujj \u3067\u5272\u3063\u3066\u3057\u307e\u3046\u306e\u3067\u3001\u30d4\u30dc\u30c3\u30c8\u9078\u629e\u306f\u5fc5\u9808…\u306a\u306e\u304b\uff1f\uff08\u3061\u3087\u3063\u3068\u3088\u304f\u308f\u304b\u3089\u305a\uff09\u3002\u6d6e\u52d5\u5c0f\u6570\u70b9\u3060\u3063\u305f\u3089\u5927\u4e08\u592b\u306a\u6c17\u3082\u3059\u308b\u306e\u3067\u3059\u304c\u3001\u3053\u308c\u306f\u3042\u3068\u3067\u78ba\u8a8d\u3057\u3066\u307f\u307e\u3059\u3002\u307e\u3042\u3001\u5272\u308a\u7b97\u306e\u3068\u3053\u308d\u306b\u30c1\u30a7\u30c3\u30af\u3092\u5165\u308c\u308b\u306e\u306f\u6570\u5024\u8a08\u7b97\u306e\u57fa\u672c\u306a\u306e\u3067\u3001\u306a\u308b\u3079\u304f\u5272\u308a\u7b97\u3067\u8aa4\u5dee\u304c\u51fa\u306a\u3044\u3088\u3046\u306b\u3059\u308b\u65b9\u5411\u3067\u3088\u3044\u304b\u3068\u3002\n<\/p>\n

\n\u3061\u306a\u307f\u306b F# MathProvider \u3092\u4f7f\u3046\u3068\n<\/p>\n

[code lang=”csharp”]
\nlet (p,L,U) = MathProvider.LinearAlgebra.lu A
\n[\/code]<\/p>\n

\n\u3068\u3057\u3066\u3001\u6b21\u306a\u308b\u7d50\u679c\u3092\u5f97\u3089\u308c\u307e\u3059\u3002\n<\/p>\n

[code lang=”csharp”]
\nval p : (int -> int)
\nval U : matrix = matrix [[2.0; 1.0; 0.0]
\n [0.0; 1.5; 1.0]
\n [0.0; 0.0; 1.666666667]]
\nval L : matrix = matrix [[1.0; 0.0; 0.0]
\n [0.5; 1.0; 0.0]
\n [0.5; 0.3333333333; 1.0]]
\n[\/code]<\/p>\n

\nMathProvider \u306e LU \u5206\u89e3\u306e\u30b3\u30fc\u30c9\u306f\u3082\u3063\u3068\u30b7\u30f3\u30d7\u30eb\u3067\u3001\u3053\u3093\u306a\u611f\u3058\u306b\u306a\u3063\u3066\u3044\u307e\u3059\u3002\u3046\u307e\u304f\u4ea4\u4e92\u306b\u8a08\u7b97\u3059\u308b\u3068 L \u306e\u307b\u3046\u306f\u307b\u3068\u3093\u3069\u8a08\u7b97\u305b\u305a\u306b\u6e08\u307f\u307f\u305f\u3044\u3067\u3059\u3002\u884c\u5217\u6570\u304c\u591a\u304f\u306a\u308b\u3068 swapR \u306e\u30d1\u30d5\u30a9\u30fc\u30de\u30f3\u30b9\u304c\u554f\u984c\u306b\u306a\u308a\u307e\u3059\u304c\u3001\u6709\u9650\u8981\u7d20\u6cd5\u4ee5\u5916\u3060\u3063\u305f\u3089\u5927\u4e08\u592b\u304b\u3068\u3002\n<\/p>\n

[code lang=”csharp”]
\nfor i=0 to nA-2 do
\n let mutable maxi = i \/\/ Find the biggest pivot element.
\n for k=i+1 to nA-1 do
\n if abs(U.[maxi,i]) < abs(U.[k,i]) then maxi <- k\n done\n \/\/let maxi = maxIndex (i+1) (nA-1) (fun k -> abs(U.[k,i]))<\/p>\n

if maxi <> i then
\n swapR U i maxi \/\/ Swap rows if needed.
\n swapR L i maxi \/\/ REVIEW can be made more performant.
\n let t = P.[i]
\n P.[i] <- P.[maxi]\n P.[maxi] <- t\n \n for j=i+1 to nA-1 do\n L.[j,i] <- U.[j,i] \/ U.[i,i]\n for k=i+1 to nA-1 do\n U.[j,k] <- U.[j,k] - L.[j,i] * U.[i,k]\n done\n U.[j,i] <- 0.0\n done\ndone\n[\/code]\n\n\n

\nsum \u306e\u3068\u3053\u308d\u306f\u3001\u3082\u3063\u3068 F# \u3089\u3057\u304f\u66f8\u304d\u63db\u3048\u3089\u308c\u306a\u3044\u304b\u306a\u3001\u3068\u601d\u6848\u4e2d\u3067\u3059\u3002\u30b3\u30fc\u30c9\u304c\u77ed\u304f\u306a\u308b\u3057\u03a3\u3092\u4f7f\u3063\u305f\u307b\u3046\u304c\u3001\u6570\u5f0f\u3092\u540c\u3058\u306b\u306a\u308b\u306e\u3067\u3002<\/p>\n","protected":false},"excerpt":{"rendered":"

\u65e7\u30bf\u30a4\u30d7\u306a\u306e\u3067\u3001\u9023\u7acb\u4e00\u6b21\u65b9\u7a0b\u5f0f\u3092\u89e3\u304f\u3068\u304d\u306b\u9006\u884c\u5217\u3092\u4f7f\u3063\u3066\u89e3\u3051\u3070OK…\u3068\u982d\u304b\u3089\u601d\u3044\u8fbc\u3093\u3067\u3044\u305f\u306e\u3067\u3059\u304c\u3001Gauss\u306e\u6d88\u53bb\u6cd5\u3092\u4f7f\u3046\u3068\u9006\u884c\u5217\u306f\u51fa\u3066\u3053\u306a\u304f\u3066\u300c\uff1f\u300d\u3068\u306a\u308b\u308f\u3051\u3067\u3059\u3002\u4e0a\u4e09\u89d2\u5316\u3092\u4f7f\u3046\u306e\u3060\u304b\u3089\u3001LU\u5206\u89e3\u3067\u3044\u3044\u3093 … \u7d9a\u304d\u3092\u8aad\u3080 →<\/span><\/a><\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"_jetpack_newsletter_access":"","_jetpack_dont_email_post_to_subs":false,"_jetpack_newsletter_tier_id":0,"_jetpack_memberships_contains_paywalled_content":false,"_jetpack_hide_featured_image":false,"_jetpack_featured_image_duplicate_dismissed":0,"_jetpack_feature_clip_id":0,"_jetpack_memberships_contains_paid_content":false,"footnotes":"","jetpack_post_was_ever_published":false},"categories":[59],"tags":[],"class_list":["post-5300","post","type-post","status-publish","format-standard","hentry","category-f"],"jetpack_sharing_enabled":true,"jetpack_featured_media_url":"","_links":{"self":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts\/5300","targetHints":{"allow":["GET"]}}],"collection":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/comments?post=5300"}],"version-history":[{"count":1,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts\/5300\/revisions"}],"predecessor-version":[{"id":5301,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/posts\/5300\/revisions\/5301"}],"wp:attachment":[{"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/media?parent=5300"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/categories?post=5300"},{"taxonomy":"post_tag","embeddable":true,"href":"http:\/\/www.moonmile.net\/blog\/wp-json\/wp\/v2\/tags?post=5300"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}