I have started this code a month ago and after struggling with optimx and the function design finally I've managed to have it up and running, but the new issue is time needed to complete the processing as it takes almost 45 mins to complete 18 parameters on DualCore 3Ghz processor. My concerns are:
1- Is this time normal for R? if so what are the other options? I was talking with a professor from the statistics and he stated I might switch to Matlab but I am not sure and I don't want to start re-writing all the code in Matlab before I make sure from the experienced users that Matlab has more computational capacity and reliability in non-linear optimization problems for the same hardware setting.
2- If this time is not normal how can I speed it up and how difficult to run the code in parallel giving that I have 0 hands on parallelization, the code is a nonlinear optimization for 18 parameters and I expect them to be more than that.
The Code:
INJ.1<-"I01 I02 I03 I04 I05
2.78E+02 1.82E+03 3.62E+02 2.90E+02 7.73E+02
7.92E+02 1.21E+03 9.33E+02 6.32E+02 5.10E+02
2.30E+03 7.54E+02 9.60E+02 6.29E+02 1.05E+03
3.61E+03 3.05E+02 7.77E+02 5.87E+02 1.02E+03
3.89E+02 1.35E+03 7.66E+02 4.00E+02 7.43E+02
1.31E+03 1.63E+03 8.95E+02 3.85E+02 1.10E+02
1.39E+03 1.16E+03 9.07E+02 4.99E+02 2.48E+02
1.94E+03 1.09E+03 8.34E+02 5.22E+02 2.48E+02
2.04E+03 1.11E+03 7.85E+02 2.67E+02 4.27E+02
1.06E+03 1.36E+03 8.80E+02 6.13E+02 7.16E+02
1.40E+03 1.29E+03 8.65E+02 6.17E+02 9.79E+02
1.20E+03 1.68E+03 6.78E+02 6.10E+02 9.30E+02
1.45E+03 1.49E+03 7.66E+02 3.81E+02 1.07E+03
1.16E+03 1.58E+03 1.09E+03 5.33E+02 8.38E+02
1.33E+03 1.38E+03 9.10E+02 6.29E+02 8.80E+02
1.33E+03 1.38E+03 9.10E+02 6.29E+02 8.80E+02
1.47E+03 1.33E+03 9.56E+02 5.79E+02 1.21E+03
1.87E+03 1.60E+03 1.08E+03 5.33E+02 1.22E+03
1.57E+03 1.76E+03 7.89E+02 5.87E+02 9.79E+02
9.68E+02 1.89E+03 8.08E+02 5.14E+02 1.08E+03
6.13E+02 1.79E+03 8.38E+02 5.37E+02 1.12E+03
1.38E+03 1.66E+03 8.19E+02 5.64E+02 9.90E+02
1.42E+03 1.55E+03 8.27E+02 5.79E+02 8.08E+02
1.35E+03 1.70E+03 9.37E+02 6.40E+02 6.82E+02
1.49E+03 1.70E+03 7.47E+02 6.86E+02 1.05E+03
7.62E+02 1.66E+03 7.01E+02 6.13E+02 7.92E+02
1.30E+03 1.62E+03 7.09E+02 6.13E+02 6.67E+02
1.47E+03 1.31E+03 6.55E+02 5.68E+02 7.16E+02
1.38E+03 1.56E+03 5.60E+02 6.97E+02 7.70E+02
1.40E+03 1.61E+03 6.82E+02 6.90E+02 7.89E+02
1.31E+03 1.58E+03 8.61E+02 7.20E+02 1.42E+03
8.80E+02 1.58E+03 8.38E+02 7.09E+02 1.46E+03
1.40E+03 1.57E+03 8.57E+02 6.74E+02 1.33E+03
1.15E+03 1.54E+03 7.85E+02 6.82E+02 1.44E+03
1.51E+03 1.53E+03 8.23E+02 6.74E+02 1.23E+03
1.08E+03 1.54E+03 8.04E+02 6.74E+02 1.19E+03
1.57E+03 1.54E+03 7.28E+02 7.01E+02 1.15E+03
1.78E+03 1.49E+03 7.85E+02 7.43E+02 1.20E+03
1.87E+03 1.46E+03 8.69E+02 8.57E+02 8.72E+02
2.14E+03 1.41E+03 1.01E+03 8.27E+02 1.78E+03
2.01E+03 1.46E+03 1.02E+03 7.73E+02 1.55E+03
1.99E+03 1.40E+03 1.11E+03 7.16E+02 1.39E+03
2.22E+03 1.44E+03 1.14E+03 7.01E+02 1.19E+03
2.11E+03 1.48E+03 9.03E+02 7.81E+02 1.21E+03
2.48E+03 1.28E+03 7.35E+02 7.28E+02 1.26E+03
2.40E+03 1.26E+03 1.07E+03 8.61E+02 1.07E+03
2.40E+03 1.17E+03 1.20E+03 8.50E+02 1.22E+03
2.43E+03 7.85E+02 1.26E+03 8.27E+02 1.07E+03
2.64E+03 6.55E+02 1.23E+03 8.15E+02 1.03E+03
2.60E+03 6.13E+02 8.27E+02 7.24E+02 9.41E+02
2.35E+03 6.74E+02 9.37E+02 6.55E+02 9.30E+02
2.77E+03 5.45E+02 1.01E+03 7.43E+02 1.14E+03
2.96E+03 6.13E+02 9.30E+02 8.42E+02 9.98E+02
3.15E+03 5.45E+02 8.95E+02 8.11E+02 1.21E+03
7.28E+02 6.06E+02 8.91E+02 8.57E+02 1.18E+03
1.59E+03 6.78E+02 8.76E+02 8.30E+02 1.23E+03
1.79E+03 4.23E+02 8.88E+02 8.19E+02 1.26E+03
4.50E+02 1.96E+03 8.38E+01 1.03E+03 1.94E+02
1.17E+03 8.65E+02 4.04E+02 7.31E+02 8.50E+02
1.68E+03 1.28E+03 1.49E+02 9.90E+02 1.58E+03
1.65E+03 1.43E+03 7.20E+02 1.19E+03 2.22E+03
1.71E+03 1.39E+03 8.30E+02 8.46E+02 2.01E+03
1.64E+03 5.83E+02 6.86E+02 9.75E+02 1.45E+03
1.57E+03 8.42E+02 5.18E+02 1.10E+03 1.19E+03
1.87E+03 1.10E+03 1.83E+02 1.13E+03 9.94E+02
1.82E+03 9.52E+02 2.32E+02 1.16E+03 1.22E+03
1.51E+03 0.00E+00 1.26E+02 1.68E+03 1.28E+03
1.39E+03 0.00E+00 2.48E+02 1.31E+03 1.10E+03
1.45E+03 0.00E+00 2.90E+02 1.24E+03 1.09E+03
1.43E+03 0.00E+00 4.19E+02 1.61E+03 1.28E+03
1.38E+03 0.00E+00 2.10E+02 1.68E+03 1.23E+03
1.07E+03 0.00E+00 5.03E+02 1.38E+03 1.21E+03
1.37E+03 0.00E+00 6.59E+02 1.44E+03 1.04E+03
1.37E+03 0.00E+00 9.45E+02 1.83E+03 1.28E+03
1.33E+03 0.00E+00 9.03E+02 1.66E+03 8.19E+02
"
INJ<-as.matrix(read.table(text=INJ.1, header=T))
PRD.1<-"P01
981.32019
1062.5702
1439.7673
1694.0723
1085.1016
1243.6089
1191.5941
1302.2167
1333.5266
1242.0212
1342.6954
1371.2767
1394.1171
1400.7926
1373.1791
1194.7649
1289.6071
1413.0968
1313.444
1243.9548
1187.6825
1225.5112
1186.6755
1195.6844
1290.3555
1107.2397
1100.14
1075.3673
1121.7092
1168.071
1377.4421
1354.8871
1375.5179
1363.9353
1355.4352
1287.1311
1317.6824
1383.2178
1365.5146
1658.0327
1625.3192
1568.1354
1542.4583
1511.9888
1493.4116
1529.6499
1590.22
1516.1698
1489.6599
1344.8517
1294.4697
1416.6488
1446.6407
1513.5117
1218.4484
1266.1337
1281.0219
915.39948
966.0614
1309.4725
1736.4828
1685.2651
1416.7249
1316.5566
1260.5198
1309.8859
1282.8658
1124.2899
1090.3422
1259.1481
1243.2527
1170.8752
1197.361
1434.2209
1291.3245
"
PRD<-as.matrix(read.table(text=PRD.1, header=T))
PR.1<-"P01 P02 P03 P04
998.5638 810 528 689
999.5638 857 985 880
1000.5638 841 552 625
1001.5638 835 742 546
1002.5638 565 745 872
1003.5638 752 957 964
1004.5638 726 987 737
1005.5638 861 695 845
1006.5638 919 938 577
1007.5638 721 631 854
1008.5638 820 702 631
1009.5638 879 648 744
1010.5638 698 820 783
1011.5638 772 679 536
1012.5638 767 785 803
1012.56378 767 785 803
1013.56378 910 775 513
1014.56378 540 828 770
1015.56378 581 723 695
1016.56378 759 580 895
1017.56378 640 863 936
1018.56378 554 735 979
1019.56378 920 748 684
1020.56378 812 613 511
1021.56378 877 851 988
1022.56378 525 894 585
1023.56378 772 534 839
1024.56378 907 631 976
1025.56378 699 724 774
1026.56378 844 647 881
1027.56378 953 885 811
1028.56378 690 776 836
1029.56378 712 660 617
1030.56378 648 669 966
1031.56378 880 525 507
1032.56378 773 700 702
1033.56378 504 541 570
1034.56378 793 860 652
1035.56378 809 766 649
1036.56378 860 847 677
1037.56378 757 909 932
1038.56378 900 594 509
1039.56378 962 987 994
1040.56378 648 958 675
1041.56378 952 971 800
1042.56378 917 900 840
1043.56378 791 553 721
1044.56378 782 615 640
1045.56378 674 563 504
1046.56378 637 988 595
1047.56378 671 809 812
1048.56378 639 786 558
1049.56378 837 610 505
1050.56378 752 950 524
1051.56378 831 807 622
1052.56378 547 949 864
1053.56378 755 839 534
1054.56378 500 517 724
1055.56378 644 507 867
1056.56378 692 804 891
1057.56378 713 726 800
1058.56378 898 952 777
1059.56378 761 587 670
1060.56378 959 847 949
1061.56378 562 667 678
1062.56378 643 863 513
1063.56378 672 996 638
1064.56378 707 890 723
1065.56378 878 794 629
1066.56378 719 512 859
1067.56378 672 944 794
1068.56378 727 963 613
1069.56378 589 831 675
1070.56378 621 945 953
1071.56378 561 796 563
"
PR<-as.matrix(read.table(text=PR.1, header=T))
lambda=c(0.0865382447763093, 0.0548765126242383,
0.0931262173744562, 0.106388796980127,
0.14140614108626,0.6326083, 0.3616249,
0.4217125, 0.3821551, 0.5059075,1,1,1,1,1,
0.0171603825377974, 0.00298285562229438,
0.000573790646237408)
#const. i.dash
fn1 <- function (lambda){
#const. i.dash
i.dash=matrix(ncol=ncol(INJ), nrow=(nrow(INJ)))
for (i in 1:ncol(INJ)){
for (j in 1:nrow (INJ)){
temp=0
for (k in 1:j){
temp=(1/lambda[ncol(INJ)+i])*exp((k-j)/lambda[ncol(INJ)+i])*INJ[k,i]+temp
}
i.dash[j,i]=temp
}
}
### preparing the i dash X lambda matrix
lambda.i=matrix(ncol=ncol(INJ),nrow=nrow(INJ))
for (i in 1: ncol(INJ)){
lambda.i[,i]=lambda[i+1]*i.dash[,i]
}
q.hat=matrix(1,nrow=nrow(i.dash), data=rowSums(lambda.i))
#########################constructing BHP matrix #####################
#######################preparing the BHP matrix#######################
##Preparing the 1st term (Pwf(n0)e)
bhp.1term=matrix(ncol=ncol(PR), nrow=nrow(PR))
for ( k in 1: ncol(PR))
for (i in 1: nrow(PR)){
bhp.1term[i,k]=PR[1,k]*exp((1-i)/lambda[((2*ncol(INJ))+k)])
}
##Preparing the 2nd term (Pwfkj)
bhp.2term=matrix(ncol=ncol(PR), nrow=nrow(PR))
bhp.2term=PR
bhp.3term=matrix(ncol=ncol(PR), nrow=nrow(PR))
## Preparing the 3rd term(P'wfkj)
for ( k in 1: ncol(PR)){
for (l in 1: nrow(PR)){
temp=0
for (m in 1:l){
temp=(1/lambda[((2*ncol(INJ))+k)]*exp((m-l)/lambda[((2*ncol(INJ))+k)])*PR[m,k]+temp)
}
bhp.3term[l,k]=temp
}
}
#################### wrapping up ##############################
bhp.master=matrix(ncol=ncol(PR), nrow=nrow(PR))
bhp.master=bhp.3term+bhp.2term+bhp.1term
BHP.nue=matrix(ncol=ncol(PR), nrow=nrow(PR), data =0)
for ( i in 1: nrow (PR)){
for ( j in 1: ncol (PR)){
BHP.nue[i,j]=bhp.master[i,j]*lambda[(2*ncol(INJ))+ncol(PR)+j]+BHP.nue[i,j]
}
}
q.hat=q.hat+rowSums(BHP.nue)
#################################################################
mini=(PRD-q.hat)^2
final=sum(mini)
return(final)
}
ini=lambda
mm=fn1(ini)
lamb=optimx(ini,fn1,hessian=T,itnmax=100,method=c("BFGS"), control=list(trace=5))
I am having a similar problem with the speed of optim() and it is taking much longer than 45 minutes - which doesn't seem too bad but I will try to help you out anyway!
I don't believe that you can an optimization function in parallel in any case since it takes steps sequentially towards a supposed max or min. If you figure that out, let me know!
I just saw this post about using the different optimization methods with optimx() and this might speed the process up. Have you tried this:
https://stat.ethz.ch/pipermail/r-help/2011-February/267999.html
If that doesn't help, I would suggest reading the Performance Criteria and Nonconvergence section of this IBM article:
http://www.ibm.com/developerworks/library/ba-optimR-john-nash/
It suggests this:
"You want to spend your effort on speeding up the objective function. It is sometimes worth putting the objective function in Fortran. The optim() function can accept a .Call expression to compiled code instead of an R function (to evaluate the objective function). Write the optimizer in R so everyone can see the implementation. The R package that is called microbenchmark runs something like a hundred times and gives you the distribution of times taken. This information is useful if you want to see where a problem is taking time."
R can compile code in Fortran and C which will run much faster. This is a nice trick for speeding up FOR loops in R and they are saying that writing the function you are optimizing in either of these languages will likely speed it up.