R version 2.9.0 (2009-04-17) Copyright (C) 2009 The R Foundation for Statistical Computing ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. Natural language support but running in an English locale R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > # orthit.r > # > # orthogonal iteration for largest eigenvalues & vectors > # > # uses modified Gram-Schmidt > mgs <- function(X) { + R <- matrix(0,ncol(X),ncol(X)) # matrices of zeros to start + Q <- matrix(0,nrow(X),ncol(X)) # + for (k in 1:ncol(X)) { + s <- sqrt( sum(X[,k]^2) ) # norm + R[k,k] <- s + Q[,k] <- X[,k]/s # normalize + if(k == ncol(X)) next + l <-((k+1):ncol(X)) + R[k,l] <- t(Q[,k])%*%X[,l] # get coefficients + X[,l] <- X[,l]-outer(Q[,k],R[k,l])} # residuals + mgs <- list(Q,R)} # return a list > # > X <- cbind(rep(1,6),(1:6),(1:6)^2) > X [,1] [,2] [,3] [1,] 1 1 1 [2,] 1 2 4 [3,] 1 3 9 [4,] 1 4 16 [5,] 1 5 25 [6,] 1 6 36 > mgsX <- mgs(X) > mgsX[[1]] [,1] [,2] [,3] [1,] 0.4082483 -0.5976143 0.5455447 [2,] 0.4082483 -0.3585686 -0.1091089 [3,] 0.4082483 -0.1195229 -0.4364358 [4,] 0.4082483 0.1195229 -0.4364358 [5,] 0.4082483 0.3585686 -0.1091089 [6,] 0.4082483 0.5976143 0.5455447 > mgsX[[2]] [,1] [,2] [,3] [1,] 2.449490 8.573214 37.150594 [2,] 0.000000 4.183300 29.283101 [3,] 0.000000 0.000000 6.110101 > # assign matrix > print("first matrix") # sparse transition matrix [1] "first matrix" > A <- matrix( c(.9,0,.1,0,0,0,0, .3,.6,.0,.1,0,0,0, + 0.,.2,.7,0,.1,0,0, 0,0,.3,.6,0,.1,0, 0,0,0,.2,.7,0,.1, + 0,0,0,0,.3,.7,0, 0,0,0,0,0,.2,.8) ,7,7) > A # print it out [,1] [,2] [,3] [,4] [,5] [,6] [,7] [1,] 0.9 0.3 0.0 0.0 0.0 0.0 0.0 [2,] 0.0 0.6 0.2 0.0 0.0 0.0 0.0 [3,] 0.1 0.0 0.7 0.3 0.0 0.0 0.0 [4,] 0.0 0.1 0.0 0.6 0.2 0.0 0.0 [5,] 0.0 0.0 0.1 0.0 0.7 0.3 0.0 [6,] 0.0 0.0 0.0 0.1 0.0 0.7 0.2 [7,] 0.0 0.0 0.0 0.0 0.1 0.0 0.8 > Z <- (diag(1,7))[,1:3] # starting orthonormal vectors > for (i in (1:40)) { + V <- A %*% Z # multiply + qrnew <- mgs(V) # normalize via MGS + Z <- qrnew[[1]] + if( (i %% 5) == 0 ) { # print a few + print(i) + print(Z) + print(qrnew[[2]])} } [1] 5 [,1] [,2] [,3] [1,] 0.934971745 -0.15599219 -0.27510557 [2,] 0.120778186 0.68766340 0.04756984 [3,] 0.324844985 0.12701798 0.65878137 [4,] 0.030463075 0.64200840 -0.18585320 [5,] 0.068534460 0.01893606 0.60588408 [6,] 0.003252179 0.27128677 -0.16944508 [7,] 0.008980789 -0.02284841 0.24023442 [,1] [,2] [,3] [1,] 0.9434068 0.2769921 0.005167419 [2,] 0.0000000 0.6434897 0.257612924 [3,] 0.0000000 0.0000000 0.677810245 [1] 10 [,1] [,2] [,3] [1,] 0.87813884 -0.39142308 0.02339387 [2,] 0.19471930 0.28672338 -0.05109485 [3,] 0.39751630 0.38151437 -0.11946461 [4,] 0.10026558 0.45748316 -0.09017657 [5,] 0.14065587 0.52322969 0.02331221 [6,] 0.03671665 0.35852824 0.03761166 [7,] 0.04185681 0.08614991 0.98614418 [,1] [,2] [,3] [1,] 0.9737193 0.1241984 -0.01010691 [2,] 0.0000000 0.8823304 0.09645843 [3,] 0.0000000 0.0000000 0.79110307 [1] 15 [,1] [,2] [,3] [1,] 0.83661416 -0.4463913 0.20914498 [2,] 0.21602836 0.1666566 -0.26698258 [3,] 0.42795096 0.2792857 -0.39280085 [4,] 0.14177261 0.4082131 -0.28095132 [5,] 0.19675604 0.5710513 0.02821056 [6,] 0.07693209 0.3618391 0.47575128 [7,] 0.07440107 0.2670109 0.65162674 [,1] [,2] [,3] [1,] 0.9846345 0.0818811 -0.02627627 [2,] 0.0000000 0.9124594 0.13635195 [3,] 0.0000000 0.0000000 0.82971812 [1] 20 [,1] [,2] [,3] [1,] 0.8030286 -0.4791282 0.3084791 [2,] 0.2241438 0.1045694 -0.4372804 [3,] 0.4430996 0.1975968 -0.5126211 [4,] 0.1691118 0.3600719 -0.2695841 [5,] 0.2402526 0.5587833 0.1423670 [6,] 0.1092587 0.4046293 0.3973374 [7,] 0.1015371 0.3388828 0.4472416 [,1] [,2] [,3] [1,] 0.9902525 0.06475123 -0.05373603 [2,] 0.0000000 0.92611017 0.12036442 [3,] 0.0000000 0.00000000 0.79377556 [1] 25 [,1] [,2] [,3] [1,] 0.7758138 -0.50526955 0.3376591 [2,] 0.2278846 0.06563388 -0.5224195 [3,] 0.4515506 0.14247129 -0.5019667 [4,] 0.1881616 0.33286615 -0.2657835 [5,] 0.2721265 0.55116335 0.1497431 [6,] 0.1335839 0.41997572 0.3333956 [7,] 0.1223979 0.35935539 0.3960771 [,1] [,2] [,3] [1,] 0.993651 0.05370988 -0.0700701 [2,] 0.000000 0.93073976 0.1124810 [3,] 0.000000 0.00000000 0.7657436 [1] 30 [,1] [,2] [,3] [1,] 0.7545155 -0.52908918 0.3456879 [2,] 0.2298167 0.04405273 -0.5406248 [3,] 0.4565102 0.11177327 -0.4883029 [4,] 0.2015333 0.31832664 -0.2794734 [5,] 0.2949836 0.54255889 0.1404568 [6,] 0.1513174 0.42179627 0.3089510 [7,] 0.1376992 0.36334095 0.3956017 [,1] [,2] [,3] [1,] 0.9958106 0.04642391 -0.07373393 [2,] 0.0000000 0.93176932 0.10992554 [3,] 0.0000000 0.00000000 0.75803252 [1] 35 [,1] [,2] [,3] [1,] 0.7384072 -0.54853171 0.3487172 [2,] 0.2308743 0.03258788 -0.5366012 [3,] 0.4595120 0.09386204 -0.4861757 [4,] 0.2109083 0.30974786 -0.2870047 [5,] 0.3111681 0.53461415 0.1304782 [6,] 0.1639949 0.41969764 0.3037796 [7,] 0.1486555 0.36240692 0.4030016 [,1] [,2] [,3] [1,] 0.997207 0.04168006 -0.07120315 [2,] 0.000000 0.93158370 0.10907645 [3,] 0.000000 0.00000000 0.75855314 [1] 40 [,1] [,2] [,3] [1,] 0.7265362 -0.56313533 0.3501848 [2,] 0.2314777 0.02602112 -0.5314705 [3,] 0.4613770 0.08253902 -0.4874900 [4,] 0.2174612 0.30429779 -0.2896718 [5,] 0.3225378 0.52823104 0.1249058 [6,] 0.1729544 0.41701670 0.3049079 [7,] 0.1564010 0.36036191 0.4059351 [,1] [,2] [,3] [1,] 0.9981211 0.03848999 -0.06801597 [2,] 0.0000000 0.93113231 0.10905378 [3,] 0.0000000 0.00000000 0.75992386 > print("second matrix") [1] "second matrix" > A <- matrix(c(1,-3,-2,1,-3,10,-3,6,-2,-3,3,-2,1,6,-2,1),4,4) > A # print it out [,1] [,2] [,3] [,4] [1,] 1 -3 -2 1 [2,] -3 10 -3 6 [3,] -2 -3 3 -2 [4,] 1 6 -2 1 > Z <- (diag(1,4))[,1:3] # starting orthonormal vectors > for (i in (1:40)) { + V <- A %*% Z # multiply + qrnew <- mgs(V) # normalize via MGS + Z <- qrnew[[1]] + if( (i %% 5) == 0 ) { # print a few + print(i) + print(Z) + print(qrnew[[2]])} } [1] 5 [,1] [,2] [,3] [1,] 0.1278801 0.4950670 -0.694874355 [2,] -0.8611733 -0.3241005 -0.391574214 [3,] 0.2676858 -0.7199610 0.002287387 [4,] -0.4127608 0.3626619 0.603169987 [,1] [,2] [,3] [1,] 14.31096 -0.6934425 0.02377826 [2,] 0.00000 4.3669807 -0.57159748 [3,] 0.00000 0.0000000 3.48997338 [1] 10 [,1] [,2] [,3] [1,] 0.1193779 0.6915587 0.5017766 [2,] -0.8559982 -0.1968821 0.4780121 [3,] 0.2794007 -0.6793546 0.2245317 [4,] -0.4182718 0.1464967 -0.6850621 [,1] [,2] [,3] [1,] 14.32951 -0.00194321 1.922750e-05 [2,] 0.00000 4.44965602 -1.634954e-01 [3,] 0.00000 0.00000000 3.420695e+00 [1] 15 [,1] [,2] [,3] [1,] 0.1193460 0.6456311 -0.5596440 [2,] -0.8559888 -0.2374806 -0.4592137 [3,] 0.2794324 -0.6961989 -0.1650047 [4,] -0.4182789 0.2051112 0.6698469 [,1] [,2] [,3] [1,] 14.32951 -5.640349e-06 1.483835e-08 [2,] 0.00000 4.456446e+00 -4.343872e-02 [3,] 0.00000 0.000000e+00 3.415483e+00 [1] 20 [,1] [,2] [,3] [1,] 0.1193459 0.6582618 0.5447313 [2,] -0.8559888 -0.2269159 0.4645250 [3,] 0.2794325 -0.6922429 0.1808844 [4,] -0.4182789 0.1897374 -0.6743628 [,1] [,2] [,3] [1,] 14.32951 -1.641555e-08 1.141265e-11 [2,] 0.00000 4.456923e+00 -1.147819e-02 [3,] 0.00000 0.000000e+00 3.415118e+00 [1] 25 [,1] [,2] [,3] [1,] 0.1193459 0.6549586 -0.5486985 [2,] -0.8559888 -0.2297184 -0.4631455 [3,] 0.2794325 -0.6933232 -0.1766986 [4,] -0.4182789 0.1938083 0.6732041 [,1] [,2] [,3] [1,] 14.32951 -4.778516e-11 8.437695e-15 [2,] 0.00000 4.456957e+00 -3.031821e-03 [3,] 0.00000 0.000000e+00 3.415092e+00 [1] 30 [,1] [,2] [,3] [1,] 0.1193459 0.6558334 0.5476526 [2,] -0.8559888 -0.2289789 0.4635115 [3,] 0.2794325 -0.6930403 0.1778048 [4,] -0.4182789 0.1927337 -0.6735125 [,1] [,2] [,3] [1,] 14.32951 -1.391665e-13 -1.110223e-16 [2,] 0.00000 4.456959e+00 -8.007965e-04 [3,] 0.00000 0.000000e+00 3.415090e+00 [1] 35 [,1] [,2] [,3] [1,] 0.1193459 0.6556025 -0.5479290 [2,] -0.8559888 -0.2291743 -0.4634150 [3,] 0.2794325 -0.6931152 -0.1775127 [4,] -0.4182789 0.1930176 0.6734312 [,1] [,2] [,3] [1,] 14.32951 -1.665335e-16 -7.771561e-16 [2,] 0.00000 4.456959e+00 -2.115144e-04 [3,] 0.00000 0.000000e+00 3.415090e+00 [1] 40 [,1] [,2] [,3] [1,] 0.1193459 0.6556635 0.5478560 [2,] -0.8559888 -0.2291227 0.4634405 [3,] 0.2794325 -0.6930954 0.1775899 [4,] -0.4182789 0.1929426 -0.6734527 [,1] [,2] [,3] [1,] 14.32951 -3.885781e-16 -2.220446e-16 [2,] 0.00000 4.456959e+00 -5.586730e-05 [3,] 0.00000 0.000000e+00 3.415090e+00 > print("third matrix") [1] "third matrix" > A <- matrix(c(8,-3,-2,0,6,5, -3,6,1,-6,0,-2, -2,1,5,-5,2,0, + 0,-6,-5,8,-3,-2, 6,0,2,-3,6,1, 5,-2,0,-2,1,5),6,6) > A # print it out [,1] [,2] [,3] [,4] [,5] [,6] [1,] 8 -3 -2 0 6 5 [2,] -3 6 1 -6 0 -2 [3,] -2 1 5 -5 2 0 [4,] 0 -6 -5 8 -3 -2 [5,] 6 0 2 -3 6 1 [6,] 5 -2 0 -2 1 5 > Z <- (diag(1,6))[,1:3] # starting orthonormal vectors > for (i in (1:40)) { + V <- A %*% Z # multiply + qrnew <- mgs(V) # normalize via MGS + Z <- qrnew[[1]] + if( (i %% 5) == 0 ) { # print a few + print(i) + print(Z) + print(qrnew[[2]])} } [1] 5 [,1] [,2] [,3] [1,] 7.219072e-01 0.0005990716 -0.23334919 [2,] -3.032295e-01 0.4549822918 -0.61624099 [3,] -7.579253e-02 0.4171729814 0.74681593 [4,] -3.438774e-17 -0.7219005897 0.06368854 [5,] 4.534964e-01 0.3036030927 0.03641670 [6,] 4.189254e-01 0.0751143097 0.05175677 [,1] [,2] [,3] [1,] 16.14272 -0.0002631331 0.02801407 [2,] 0.00000 16.1421134278 0.13667093 [3,] 0.00000 0.0000000000 4.45616341 [1] 10 [,1] [,2] [,3] [1,] 0.72178249 9.700932e-07 -0.23355617 [2,] -0.30355342 4.534231e-01 -0.61753532 [3,] -0.07541086 4.190535e-01 0.74566469 [4,] 0.00000000 -7.217825e-01 0.06512224 [5,] 0.45342074 3.035540e-01 0.03520440 [6,] 0.41905631 7.540979e-02 0.05104470 [,1] [,2] [,3] [1,] 16.14274 -6.755805e-10 0.0000449000 [2,] 0.00000 1.614274e+01 0.0002190407 [3,] 0.00000 0.000000e+00 4.4559896387 [1] 15 [,1] [,2] [,3] [1,] 7.217823e-01 1.554627e-09 -0.23355776 [2,] -3.035539e-01 4.534206e-01 -0.61753509 [3,] -7.541026e-02 4.190565e-01 0.74566386 [4,] -4.126521e-17 -7.217823e-01 0.06512748 [5,] 4.534206e-01 3.035539e-01 0.03520405 [6,] 4.190565e-01 7.541026e-02 0.05104587 [,1] [,2] [,3] [1,] 16.14274 -9.992007e-16 7.195820e-08 [2,] 0.00000 1.614274e+01 3.510417e-07 [3,] 0.00000 0.000000e+00 4.455990e+00 [1] 20 [,1] [,2] [,3] [1,] 7.217823e-01 2.491509e-12 -0.23355776 [2,] -3.035539e-01 4.534206e-01 -0.61753511 [3,] -7.541026e-02 4.190565e-01 0.74566385 [4,] -2.063261e-17 -7.217823e-01 0.06512747 [5,] 4.534206e-01 3.035539e-01 0.03520404 [6,] 4.190565e-01 7.541026e-02 0.05104585 [,1] [,2] [,3] [1,] 16.14274 -1.110223e-16 1.153227e-10 [2,] 0.00000 1.614274e+01 5.625927e-10 [3,] 0.00000 0.000000e+00 4.455990e+00 [1] 25 [,1] [,2] [,3] [1,] 0.72178228 3.972941e-15 -0.23355776 [2,] -0.30355394 4.534206e-01 -0.61753511 [3,] -0.07541026 4.190565e-01 0.74566385 [4,] 0.00000000 -7.217823e-01 0.06512747 [5,] 0.45342063 3.035539e-01 0.03520404 [6,] 0.41905653 7.541026e-02 0.05104585 [,1] [,2] [,3] [1,] 16.14274 6.661338e-16 1.849354e-13 [2,] 0.00000 1.614274e+01 9.007413e-13 [3,] 0.00000 0.000000e+00 4.455990e+00 [1] 30 [,1] [,2] [,3] [1,] 7.217823e-01 -8.402851e-18 -0.23355776 [2,] -3.035539e-01 4.534206e-01 -0.61753511 [3,] -7.541026e-02 4.190565e-01 0.74566385 [4,] 6.877536e-18 -7.217823e-01 0.06512747 [5,] 4.534206e-01 3.035539e-01 0.03520404 [6,] 4.190565e-01 7.541026e-02 0.05104585 [,1] [,2] [,3] [1,] 16.14274 1.110223e-16 2.914335e-16 [2,] 0.00000 1.614274e+01 1.186551e-15 [3,] 0.00000 0.000000e+00 4.455990e+00 [1] 35 [,1] [,2] [,3] [1,] 7.217823e-01 2.274357e-18 -0.23355776 [2,] -3.035539e-01 4.534206e-01 -0.61753511 [3,] -7.541026e-02 4.190565e-01 0.74566385 [4,] 6.877536e-18 -7.217823e-01 0.06512747 [5,] 4.534206e-01 3.035539e-01 0.03520404 [6,] 4.190565e-01 7.541026e-02 0.05104585 [,1] [,2] [,3] [1,] 16.14274 -6.661338e-16 3.608225e-16 [2,] 0.00000 1.614274e+01 -8.604228e-16 [3,] 0.00000 0.000000e+00 4.455990e+00 [1] 40 [,1] [,2] [,3] [1,] 7.217823e-01 -6.489399e-18 -0.23355776 [2,] -3.035539e-01 4.534206e-01 -0.61753511 [3,] -7.541026e-02 4.190565e-01 0.74566385 [4,] 6.877536e-18 -7.217823e-01 0.06512747 [5,] 4.534206e-01 3.035539e-01 0.03520404 [6,] 4.190565e-01 7.541026e-02 0.05104585 [,1] [,2] [,3] [1,] 16.14274 2.220446e-16 7.632783e-16 [2,] 0.00000 1.614274e+01 9.159340e-16 [3,] 0.00000 0.000000e+00 4.455990e+00 > print("fourth matrix") [1] "fourth matrix" > A <- matrix(0,5,5) > A <- dbinom(row(A)-1,4,prob=(2*col(A)-1)/10) > A # print it out [,1] [,2] [,3] [,4] [,5] [1,] 0.6561 0.2401 0.0625 0.0081 0.0001 [2,] 0.2916 0.4116 0.2500 0.0756 0.0036 [3,] 0.0486 0.2646 0.3750 0.2646 0.0486 [4,] 0.0036 0.0756 0.2500 0.4116 0.2916 [5,] 0.0001 0.0081 0.0625 0.2401 0.6561 > Z <- (diag(1,5))[,1:3] # starting orthonormal vectors > for (i in (1:40)) { + V <- A %*% Z # multiply + qrnew <- mgs(V) # normalize via MGS + Z <- qrnew[[1]] + if( (i %% 5) == 0 ) { # print a few + print(i) + print(Z) + print(qrnew[[2]])} } [1] 5 [,1] [,2] [,3] [1,] 0.6498808 -0.46450354 0.4469193 [2,] 0.5741349 -0.09465529 -0.2539287 [3,] 0.4001446 0.30970852 -0.5477300 [4,] 0.2597850 0.56304880 -0.0739913 [5,] 0.1428993 0.60194145 0.6559729 [,1] [,2] [,3] [1,] 0.9533389 0.1619047 0.03654579 [2,] 0.0000000 0.8068338 0.18308686 [3,] 0.0000000 0.0000000 0.49749750 [1] 10 [,1] [,2] [,3] [1,] 0.5064127 -0.53885695 0.5295820 [2,] 0.5169924 -0.29538207 -0.2005289 [3,] 0.4544863 0.07832866 -0.5862212 [4,] 0.4046875 0.44182655 -0.1887240 [5,] 0.3254770 0.64887329 0.5477764 [,1] [,2] [,3] [1,] 0.993719 0.06600125 0.02311734 [2,] 0.000000 0.80482702 0.01670677 [3,] 0.000000 0.00000000 0.48013639 [1] 15 [,1] [,2] [,3] [1,] 0.4497441 -0.57972938 0.5374020 [2,] 0.4843182 -0.34896005 -0.1957900 [3,] 0.4596189 0.02364056 -0.5869754 [4,] 0.4471412 0.39591594 -0.1948969 [5,] 0.3898478 0.62034863 0.5388325 [,1] [,2] [,3] [1,] 0.9993028 0.02234270 0.022912854 [2,] 0.0000000 0.80055561 0.002162020 [3,] 0.0000000 0.00000000 0.480001500 [1] 20 [,1] [,2] [,3] [1,] 0.4300758 -0.593882859 0.5380460 [2,] 0.4723488 -0.365208537 -0.1953959 [3,] 0.4601413 0.007579856 -0.5869932 [4,] 0.4601531 0.380495450 -0.1953267 [5,] 0.4104272 0.607529121 0.5381574 [,1] [,2] [,3] [1,] 0.9999248 0.007353496 0.0229270214 [2,] 0.0000000 0.800060099 0.0004540977 [3,] 0.0000000 0.000000000 0.4800000264 [1] 25 [,1] [,2] [,3] [1,] 0.4235210 -0.598529152 0.5380970 [2,] 0.4683051 -0.370396123 -0.1953644 [3,] 0.4601965 0.002470457 -0.5869937 [4,] 0.4643084 0.375396879 -0.1953590 [5,] 0.4170818 0.603026240 0.5381056 [,1] [,2] [,3] [1,] 0.999992 0.002411216 0.0229294701 [2,] 0.000000 0.800006460 0.0001290338 [3,] 0.000000 0.000000000 0.4800000006 [1] 30 [,1] [,2] [,3] [1,] 0.4213618 -0.6000476148 0.5381010 [2,] 0.4669676 -0.3720822984 -0.1953619 [3,] 0.4602024 0.0008084719 -0.5869937 [4,] 0.4656579 0.3737202631 -0.1953615 [5,] 0.4192518 0.6015231787 0.5381016 [,1] [,2] [,3] [1,] 0.9999991 0.0007902022 2.292976e-02 [2,] 0.0000000 0.8000006938 4.074448e-05 [3,] 0.0000000 0.0000000000 4.800000e-01 [1] 35 [,1] [,2] [,3] [1,] 0.4206531 -0.6005443915 0.5381013 [2,] 0.4665279 -0.3726334568 -0.1953617 [3,] 0.4602030 0.0002648382 -0.5869937 [4,] 0.4660988 0.3731701311 -0.1953617 [5,] 0.4199616 0.6010280557 0.5381013 [,1] [,2] [,3] [1,] 1 0.0002589397 2.292979e-02 [2,] 0 0.8000000745 1.323159e-05 [3,] 0 0.0000000000 4.800000e-01 [1] 40 [,1] [,2] [,3] [1,] 0.4204207 -0.6007070630 0.5381013 [2,] 0.4663837 -0.3728139225 -0.1953617 [3,] 0.4602031 0.0000867758 -0.5869937 [4,] 0.4662431 0.3729897757 -0.1953617 [5,] 0.4201941 0.6008655618 0.5381013 [,1] [,2] [,3] [1,] 1 8.484982e-05 2.292979e-02 [2,] 0 8.000000e-01 4.326430e-06 [3,] 0 0.000000e+00 4.800000e-01 > rm(list=ls()) # done > q()