r语言 - 是否有任何方法可以加快此代码的最佳处理时间



我一个月前开始了这段代码,在与optimx和函数设计作斗争之后,我终于设法启动并运行了它,但新的问题是完成处理所需的时间,因为在 DualCore 45Ghz 处理器上完成 18 个参数需要将近 38 分钟。我担心的是:

1-这个时间对R来说是正常的吗?如果是这样,还有什么其他选择?我正在与一位统计学教授交谈,他说我可能会切换到Matlab但我不确定,在我从有经验的用户那里确保Matlab在非线性优化问题中具有更多的计算能力和可靠性之前,我不想开始重写 Matlab 中的所有代码相同的硬件设置。

2-如果这个时间不正常,我该如何加快速度以及并行运行代码的难度有多大,因为我在并行化方面有0手,代码是18个参数的非线性优化,我希望它们不止于此。

守则:

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))

我在 optim() 的速度方面遇到了类似的问题,它花费的时间比 45 分钟长得多 - 这似乎还不错,但无论如何我都会尽力帮助你!

我不相信在任何情况下您都可以并行执行优化函数,因为它会按顺序朝着假定的最大值或最小值迈出步进。如果您弄清楚了,请告诉我!

我刚刚看到这篇关于使用 optimx() 使用不同优化方法的文章,这可能会加快这个过程。你试过这个吗:

https://stat.ethz.ch/pipermail/r-help/2011-February/267999.html

如果这没有帮助,我建议阅读这篇 IBM 文章的性能标准和非收敛部分:

http://www.ibm.com/developerworks/library/ba-optimR-john-nash/

它表明:

"你想把精力花在加速目标函数上。有时值得将目标函数放在Fortran中。optim() 函数可以接受 .调用表达式以编译代码而不是 R 函数(以计算目标函数)。编写优化器,以便每个人都可以看到实现。称为微基准测试的 R 包运行大约一百次,并为您提供所花费时间的分布。如果您想查看问题需要花费时间的地方,此信息非常有用。

R可以用Fortran和C编译代码,这将运行得更快。这是在R中加速FOR循环的一个很好的技巧,他们说用这两种语言中的任何一种编写您正在优化的函数都可能会加快速度。

最新更新