I haven't tried running your code, but a quick read suggests you should
1) set up the input data so your code can be run with source() without any preprocessing.
2) compute the function for several sets of parameters to make sure it is correct. Maybe
create a very simple test case you can more or less hand calculate as a check.
3) rescale the problem. You have the 1st parameter very large compared to the second
4) provide an analytic gradient -- item (3) suggests that the difference in scale of
the parameters could be causing very poor gradient approximations. It's quite possible
that optim() uses a quite simple forward difference approximation. For your problem,
I believe the analytic gradient is "simple" but very tedious.
5) once scaled, you could probably get away with using nmkb from dfoptim, which is a
Nelder-Mead method but with a form of bounds.
6) You might be able to use control$parscale to do the scaling. Do so carefully -- it
is easy to get the scaling the "wrong way round".
FYI, and that of others, there is a new package optimrx on r-forge (optimr on CRAN is
also available, but has few optimizers -- it is there so the name stays available).
http://download.r-forge.r-project.org/src/contrib/optimrx_2016-9.16.tar.gz or
http://download.r-forge.r-project.org/bin/windows/contrib/latest/optimrx_2016-9.16.zip
or https://r-forge.r-project.org/R/?group_id=395 for the project and Subversion access.
This has a function optimr() that uses the optim() syntax throughout for about 20
methods, and parscale control is available for all.
John Nash
On 16-09-29 03:53 PM, Narendra Modi wrote:
> I have put together a R snippet wherein I am trying to get optimum
> values for which error is minimized. The error is the difference
> between two matrices.
> Every time I run the below code, I don't see any optimization
> happening as in the final answer is the same as the initial estimate
> regardless of what I mention as initial estimate.
> Could you please advise on how can I improve it?
> I have also tried using nloptr but to no avail.
>
> To optimize vp.kval
>
>
>
>
> my.data.matrix.prod <- matrix(a,nrow = length(a),1)
> Estimated.Qt.mat <- matrix(b,nrow = nrow(my.data.matrix.prod), ncol = 1)
> Cum.WInj.matrix <- matrix (c,nrow = nrow(my.data.matrix.prod), ncol = 1)
> Koval.tD <- matrix(,nrow = nrow(my.data.matrix.prod), ncol = 1) # tD Matrix
> Koval.fw <- matrix(,nrow = nrow(my.data.matrix.prod), ncol = 1) # fw Matrix
>
> Koval.Error.func <- function(vp.kval,n)
> {
> #First convert vector(Koval.InitialData.list) to MATRIX and send it
> to calculate estimated matrix
>
> Koval.InitialData.Matrix <- matrix(vp.kval,nrow = 2, 1,byrow = TRUE)
> # Define Koval Parameters Matrix for the "n"
>
> Qo.Koval <- Qo.Koval(Koval.InitialData.Matrix) # Get Qo Estimation from Koval
>
> diff.values <- my.data.matrix.prod[,n] - Qo.Koval #FIND DIFFERENCE
> BETWEEN CAL. MATRIX AND ORIGINAL MATRIX
>
> Error <- ((colSums ((diff.values^2), na.rm = TRUE, dims =
> 1))/nrow(my.data.matrix.prod))^0.5 #sum of square root of the diff
>
> Error # return error value
> }
>
> Qo.Koval <- function(Koval.InitialData.Matrix)
> {
> Qo.Koval.Est <- matrix(,nrow(my.data.matrix.prod), 1)
> #ncol(my.data.matrix.prod)
>
> for(rowno in 1:nrow(my.data.matrix.prod)) #number of rows of data
> {
> Koval.tD[rowno,1] = Cum.WInj.matrix[rowno,1] *
> Koval.InitialData.Matrix[1,1] # Evaluate tD matrix
>
> if(Koval.tD[rowno,1] < (1/Koval.InitialData.Matrix[2,1]))
> {
> Koval.fw[rowno,1] = 0
> }
> else
> {
> if(Koval.tD[rowno,1] > Koval.InitialData.Matrix[2,1])
> {
> Koval.fw[rowno,1] = 1
> }else
> {
> Koval.fw[rowno,1] = (Koval.InitialData.Matrix[2,1] -
> sqrt(Koval.InitialData.Matrix[2,1]/Koval.tD[rowno,1]))/(Koval.InitialData.Matrix[2,1]-1)
> }
> }
>
> Qo.Koval.Est[rowno,1] = Koval.fw[rowno,1] * Estimated.Qt.mat[rowno,1]
> }
>
> Qo.Koval.Est # Return Estimated matrix
> }
>
> vp.kval <- c(100000,1) #initial estimate
> Koval.lb = c(0,0)
> Koval.ub = c(Inf,Inf)
> n <- 1
> Koval.result <- optim(vp.kval,Koval.Error.func, method = "L-BFGS-B",
> lower = Koval.lb,
> upper = Koval.ub, n=n)
> print(paste(Koval.result$convergence))
> print(paste(Koval.result$par))
>
>
> Here is the input data:
>
> a:
>
> structure(c(414, 40, 639, 616, 677, 598, 586, 494, 322, 351,
> 322, 213, 395, 358, 406, 384, 409, 404, 370, 376, 412, 404, 369,
> 391, 341, 350, 349, 313, 302, 196, 386, 330, 350, 323, 454, 465,
> 465, 399, 416, 396, 453, 388, 496, 379, 472, 491, 492, 503, 516,
> 454, 630, 547, 578, 312, 764, 672, 548, 611, 546, 552, 520, 486,
> 581, 559, 433, 262, 650, 615, 542, 571, 542, 529, 577, 469, 557,
> 540, 546, 519, 376, 605, 520, 435, 299, 531, 538, 475, 511, 487,
> 490, 494, 537, 482, 438, 498, 312, 476, 383, 382), .Dim = c(98L,
> 1L), .Dimnames = list(NULL, "Q2"))
>
> b:
>
> structure(c(2342.12883525675, 2595.06229039124, 2715.2774272809,
> 2742.14586849367, 2678.48814516156, 2769.17482063132, 2809.26904957691,
> 2647.26143288146, 2142.48588931211, 1986.26692938822, 2417.80180308667,
> 2539.99173834861, 2889.68696098066, 2949.03395956634, 3345.265659123,
> 3178.09552101488, 3202.97894028497, 3294.04615708455, 3273.96002181006,
> 3290.59294404149, 3074.57078080845, 2809.00966959208, 2870.20594457832,
> 2994.89960881099, 3031.51083818418, 2838.72778780229, 2779.83367818986,
> 2471.70302686638, 2277.52074079803, 2313.67080371772, 2415.57558854185,
> 2593.57170885689, 2579.65222779155, 2542.47630393453, 2610.16334633228,
> 2715.1622580481, 2680.04491562794, 2676.08878142995, 2890.5657368073,
> 2939.98447437336, 2932.41354171428, 2699.29100102243, 2748.9757584712,
> 2885.90115387751, 2841.03004973532, 3111.28842226602, 3293.09352655985,
> 3448.16679970445, 3470.58231818316, 3077.6191619663, 2892.81263635983,
> 2563.00601428125, 2410.40833201752, 2696.80369889632, 3250.95996536945,
> 3900.33363068933, 3571.89127039948, 3569.09158205254, 3718.94141619046,
> 3963.05018539626, 4317.67764180387, 4143.2306512351, 4482.33003541385,
> 4313.47162218783, 4162.58533919444, 4119.75974744111, 4080.01960112015,
> 4146.78116940934, 3848.98992961189, 3507.00912988581, 3336.3269842557,
> 3691.50683899193, 3616.0923981163, 3325.14304882807, 3471.79805853069,
> 3229.60965194249, 3106.05768279943, 3184.66721766981, 3140.79657087168,
> 3242.97205541341, 3090.78617601495, 3086.74973135927, 3317.74000570974,
> 3594.90929884806, 3716.02759860505, 3622.91307702134, 3793.8518218782,
> 3666.82966979173, 3779.4557494045, 3750.98605852729, 3333.45681985961,
> 3057.22984206785, 3395.04273620089, 3623.47886822009, 3690.34495906538,
> 3827.97665203175, 3933.61679986677, 3762.82354740958), .Dim = c(98L,
> 1L))
>
> c:
>
> structure(c(2854.17262019504, 91576.5893971961, 171680.262910889,
> 257565.867448752, 335812.78671975, 424624.296030534, 510229.898689908,
> 586994.432148103, 636896.230541501, 691311.784820203, 780382.614051205,
> 860628.109248455, 961649.745761829, 1055011.51571743, 1162113.22730208,
> 1255164.70993334, 1352077.97513698, 1457172.94644257, 1554726.68952114,
> 1657279.26732716, 1745523.14071769, 1821000.62788843, 1911979.2340704,
> 2005954.86455037, 2101129.54803795, 2182822.62676551, 2258661.15941603,
> 2325202.52364728, 2387098.71588595, 2460005.26984465, 2535846.63352407,
> 2622071.03988945, 2701584.84087477, 2776628.22472118, 2859757.87551639,
> 2944689.28669068, 3026621.7086177, 3109451.02390273, 3200463.68736646,
> 3293220.09008416, 3380941.82063061, 3456992.53009029, 3541106.87910663,
> 3635049.74483035, 3721653.58199201, 3823940.56521733, 3931974.7704047,
> 4040554.293988, 4148875.73758196, 4231424.94909108, 4306157.82023537,
> 4374820.38189491, 4442080.07977206, 4535051.28823047, 4650928.35668784,
> 4793084.92990124, 4893067.57046614, 5000047.61946087, 5120237.59556207,
> 5247211.61097682, 5392662.33159569, 5515394.91894634, 5652397.35652795,
> 5780590.26125026, 5900471.93425283, 6026783.31719041, 6147868.09782702,
> 6278602.62144284, 6388178.16489299, 6482065.35854026, 6579907.11054506,
> 6702412.42037957, 6812043.87236422, 6905604.01572694, 7007786.71329603,
> 7099980.68361161, 7189071.57372477, 7290368.20702837, 7383139.53472758,
> 7487014.64958016, 7577849.79802546, 7670318.64215094, 7780726.13033402,
> 7897750.56237753, 8016910.17077053, 8126173.95193153, 8241923.12215802,
> 8351438.92663496, 8468551.68021943, 8583900.77567578, 8670079.97000769,
> 8755816.49103472, 8872115.39013376, 8988383.34061776, 9104971.76148562,
> 9224368.08502439, 9349766.5439318, 9460826.05419725), .Dim = c(98L,
> 1L))
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>