Narendra Modi
2016-Sep-29 19:53 UTC
[R] Optimization issue - Solution converges for any initial value
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))
ProfJCNash
2016-Sep-29 20:47 UTC
[R] Optimization issue - Solution converges for any initial value
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. >