Rank: 4 based on 1963 downloads (last 30 days) and 23 files submitted
photo

John D'Errico

E-mail
Company/University
Consultant
Lat/Long
43.22999954223633, -76.93000030517578

Personal Profile:

http://www.youtube.com/watch?v=IqdtzJvliMk

Professional Interests:
Numerical analysis, mathematical modeling

 

Watch this Author's files

 

Files Posted by John View all
Updated   File Tags Downloads
(last 30 days)
Comments Rating
06 Feb 2012 fminsearchbnd, fminsearchcon Bound constrained optimization using fminsearch Author: John D'Errico optimization, fminsearch, bounds, constraints, nonlinear, inequalities 292 46
  • 4.97368
5.0 | 38 ratings
01 Feb 2012 Screenshot interparc Distance based interpolation along a general curve in space Author: John D'Errico interpolation, interpolant, space curve, arc length, arc, distance 42 0
01 Feb 2012 Screenshot arclength Compute the arclength of a general curve in any number of dimensions Author: John D'Errico spline, arc length, curve, space curve, interpolation, arc 38 0
01 Feb 2012 Screenshot distance2curve Find the closest point on a (n-dimensional) curve to any given point or set of points Author: John D'Errico spline, pchip, interpolation, closest point, distance, curve 56 2
  • 5.0
5.0 | 1 rating
26 Jan 2012 Screenshot A suite of minimal bounding objects Suite of tools to compute minimal bounding circles, rectangles, triangles, spheres, incircles, etc. Author: John D'Errico minimal, minimum, area, perimeter, polygon, rectangle 63 0
Comments and Ratings by John View all
Updated File Comments Rating
06 Feb 2012 fminsearchbnd, fminsearchcon Bound constrained optimization using fminsearch Author: John D'Errico

(I've fixed the bug about the outputfcn problem and re-posted it. The new version should appear sometime this morning.)

As far as Stefan's problem goes, anytime you throw problems with variables that vary by 20 to 30 powers of 10 at ANY numerical optimizer, expect serious things to go wrong. This is a no-no for literally ANY numerical tool. And setting a convergence tolerance (TolX) to eps is just as silly, especially when your parameters vary by that many orders of magnitude.

Essentially, you are asking that one variable, which can apparently vary somewhere between 1 and 1e30, must be computed to an absolute accuracy of roughly 2e-16.

A common solution is to normalize your variables, so that both are at least of similar orders of magnitude. If you were trying to solve a problem where one variable had bounds of 1e20 to 5e20, and a second variable was bounded in the interval 1e-6 to 1e-5, then normalizing the variables to both be of order 1 would make complete sense. But your variables each vary by many powers of 10!

One thing all novices need to learn about computing, is that when things vary by that many powers of 10, is to use logs! Let fminsearch vary the log10 of your variables, so they now have bound vectors that look like [0, -6] for the lower bounds, and [30, 0] for the upper bound. Now, inside your function, take the parameters and raise 10 to that power BEFORE they are used. Do the same with the output when it is returned. As far as fminsearch is concerned, your numbers are perfectly normal looking, but your objective sees the numbers in their proper scales.

You will find that any optimizer runs much more happily when you do this, as now the variables are very nicely normalized. All of the mathematics works more simply when you work in the log space. In fact, this is surely the natural way to work with variables that vary by many powers of 10. A nice consequence is the tolerances in TolX become relative tolerances on the variables, something that surely makes much more sense.

And DON'T use eps as a convergence tolerance. Fminsearch will never be able to give you even a reasonable shot at convergence in even 10000 function evals with that fine of a tolerance. So what happens is fminsearch will keep on iterating until it runs out of function evaluations. Use a meaningful, realistic tolerance. You were probably only setting that fine of a tolerance because of the ridiculous variable scaling anyway.

04 Feb 2012 inpaint_nans Interpolates (& extrapolates) NaN elements in a 2d array. Author: John D'Errico

For the answer to your questions, you need to dive more deeply into the equations, and how they arise. The factor of 2 you found actually appears to be correct and consistent in context.

If you prefer to make up any of a variety of arbitrary equations that make sense to you, of course you can get any result you like. The fact is, that will need to wait until you write your own code.

I'm sorry, but sitting here arguing with someone who has not bothered to understand the underlying code that a factor of 2 should or should not belong seems silly to me. For some reason you seem to be trying to reverse engineer the underlying equations from a few coefficients, and are then deciding if they make sense to you, without actually expending the effort. It is the wrong direction to go when you are trying to understand a piece of software. Start from the beginning, not the end.

03 Feb 2012 inpaint_nans Interpolates (& extrapolates) NaN elements in a 2d array. Author: John D'Errico

Willem - Personally, I think that 1.57982346157884756363579013421232435 is a slightly better value, although I am willing to negotiate on the last few digits. My point is, you are apparently trying to second guess the value predicted by less than 10% on a very uncertain value. How many angels can fit on the head of a pin anyway?

So feel perfectly free to pick whatever value you like there, when you write your own code. The value that was generated is the one consistent with the equations I specified, including the boundary conditions I chose. While those equations make some rational sense, I've won't claim that the predictions inpaint_nans makes are perfect, or that it will make all users blissfully happy. That would be impossible anyway, and I'm just too busy counting angels.

02 Feb 2012 Variable Precision Integer Arithmetic Arithmetic with integers of fully arbitrary size. Arrays and vectors of vpi numbers are supported. Author: John D'Errico

Maybe we are looking at this the wrong way.

Suppose your goal is to generate random numbers, WITHOUT replacement in the interval [0,2^200 - 1].

Yes, the odds of a collision are infinitesimal, but the point is as easily made for numbers in a smaller interval. Here I'll use 20 "digit" numbers, with each digit in base 1024.

n = 500000;
spares = 100;
B = floor(rand(n + spares,20)*1024);

Eliminate any collisions, although there essentially won't ever be any for numbers this large. Unique sorts them of course, so we now extract the number we originally wanted.

B = unique(B,'rows');
ind = randperm(n+spares);
B = B(ind(1:n),:);

Having done this, you can easily convert them into VPI numbers if you want them in a big decimal form for some reason, or simply leave then as they are. Or, suppose you wanted to use a decimal form, with bounds of [0,1e1000-1]. Here we could have used numbers with 100 base 1e10 digits. The same approach applies.

Generating a few spares is incredibly cheap here, so there is no real problem.

In fact, as long as you can factor the upper limit (+1) into a list of small primes, all of which are less than 2^53, one can always use a scheme as I have described above. The only problem is if you wanted to generate samples that lie in an interval like this:

 [0, 12323242124232456456345346]

Of course, I've chosen the number 12323242124232456456345347 carefully, to be a large prime number. This is why my randint code is slow, as it allows any general upper limit. (I can certainly speed that up with some thought though, since I chose an algorithm that may not have been optimal for speed. I'll need to leave that for the future re-write.)

02 Feb 2012 Variable Precision Integer Arithmetic Arithmetic with integers of fully arbitrary size. Arrays and vectors of vpi numbers are supported. Author: John D'Errico

While I'm happy to see that some find VPI useful, nothing forces you to use a tool that is designed to operate on numbers with hundreds or even many many thousands of decimal digits. If your goal really is nothing larger than that will fit into a 64 bit integer, then use uint64 for your operations. (Only 64 bits is amazingly small in context of the problems VPI is designed to handle.) VPI obviously cannot compete in speed with arithmetic on those numbers, as this is an issue of apples versus oranges.

One day I expect to get around to re-writing VPI, but even then I will not reasonably expect more than about 10-1 speedup in any tests I've done, and even that will force some potential compromises on my part.

Use the right tool for your problem. In this case, uint64 seems like a better choice.

Comments and Ratings on John's Files View all
Updated File Comment by Comments Rating
07 Feb 2012 (Block) tri-diagonal matrices Generate (block) tridiagonal matrices Author: John D'Errico Dany
06 Feb 2012 fminsearchbnd, fminsearchcon Bound constrained optimization using fminsearch Author: John D'Errico D'Errico, John

(I've fixed the bug about the outputfcn problem and re-posted it. The new version should appear sometime this morning.)

As far as Stefan's problem goes, anytime you throw problems with variables that vary by 20 to 30 powers of 10 at ANY numerical optimizer, expect serious things to go wrong. This is a no-no for literally ANY numerical tool. And setting a convergence tolerance (TolX) to eps is just as silly, especially when your parameters vary by that many orders of magnitude.

Essentially, you are asking that one variable, which can apparently vary somewhere between 1 and 1e30, must be computed to an absolute accuracy of roughly 2e-16.

A common solution is to normalize your variables, so that both are at least of similar orders of magnitude. If you were trying to solve a problem where one variable had bounds of 1e20 to 5e20, and a second variable was bounded in the interval 1e-6 to 1e-5, then normalizing the variables to both be of order 1 would make complete sense. But your variables each vary by many powers of 10!

One thing all novices need to learn about computing, is that when things vary by that many powers of 10, is to use logs! Let fminsearch vary the log10 of your variables, so they now have bound vectors that look like [0, -6] for the lower bounds, and [30, 0] for the upper bound. Now, inside your function, take the parameters and raise 10 to that power BEFORE they are used. Do the same with the output when it is returned. As far as fminsearch is concerned, your numbers are perfectly normal looking, but your objective sees the numbers in their proper scales.

You will find that any optimizer runs much more happily when you do this, as now the variables are very nicely normalized. All of the mathematics works more simply when you work in the log space. In fact, this is surely the natural way to work with variables that vary by many powers of 10. A nice consequence is the tolerances in TolX become relative tolerances on the variables, something that surely makes much more sense.

And DON'T use eps as a convergence tolerance. Fminsearch will never be able to give you even a reasonable shot at convergence in even 10000 function evals with that fine of a tolerance. So what happens is fminsearch will keep on iterating until it runs out of function evaluations. Use a meaningful, realistic tolerance. You were probably only setting that fine of a tolerance because of the ridiculous variable scaling anyway.

06 Feb 2012 fminsearchbnd, fminsearchcon Bound constrained optimization using fminsearch Author: John D'Errico Stefan

The following test shows strange results:
----------------------------------
G0bnd 2.22e+014 L0bnd 2.707e-006

G0 2e+014 L0 3e-006

----------------------------------

function test()
options=optimset('Display','on', ...
                   'MaxFunEvals',1e4, ...
                   'MaxIter',1e4, ...
                   'TolFun', eps, ...
                   'TolX', eps, ...
                   'FunValCheck', 'on' ...
                   );
 
 w= 150e-6;
 xData = 1e-6:1e-6:w;
 G0 = 2e14;
 L0 = 3e-6;
 yData = G0*exp(-xData./L0);
 
 start_params = [G0, L0];
 min_params = [ 1, 1e-6];
 max_params = [1e30, 1];
 
 
 plot(xData,yData,'-r');
 hold('on');
 
 
 result_params=fminsearchbnd(@fitG_error_function,start_params,min_params,max_params,options,xData,yData);
 ['G0bnd ' num2str(result_params(1),4) ' L0bnd ' num2str(result_params(2),4)]
 plot(xData, result_params(1)*exp(-xData./result_params(2)),'-g');
 
 result_params=fminsearch(@fitG_error_function,start_params,options,xData,yData);
 ['G0 ' num2str(result_params(1),4) ' L0 ' num2str(result_params(2),4)]
 plot(xData, result_params(1)*exp(-xData./result_params(2)),'-b');
  
 
 function fiterror = fitG_error_function(start_params,xData,yData)
   Fitted_Curve= start_params(1) * exp(-xData./start_params(2));
   Error_Vector=Fitted_Curve - yData;
   fiterror=sum(Error_Vector.^2);
   results.fiterror=fiterror;
 end
end

06 Feb 2012 fminsearchbnd, fminsearchcon Bound constrained optimization using fminsearch Author: John D'Errico Romain

Great function, very handy.

I have noticed a small error, though. If fminsearch is not called, the function returns output.funcount, whereas it returns output.funcCount otherwise.

04 Feb 2012 inpaint_nans Interpolates (& extrapolates) NaN elements in a 2d array. Author: John D'Errico D'Errico, John

For the answer to your questions, you need to dive more deeply into the equations, and how they arise. The factor of 2 you found actually appears to be correct and consistent in context.

If you prefer to make up any of a variety of arbitrary equations that make sense to you, of course you can get any result you like. The fact is, that will need to wait until you write your own code.

I'm sorry, but sitting here arguing with someone who has not bothered to understand the underlying code that a factor of 2 should or should not belong seems silly to me. For some reason you seem to be trying to reverse engineer the underlying equations from a few coefficients, and are then deciding if they make sense to you, without actually expending the effort. It is the wrong direction to go when you are trying to understand a piece of software. Start from the beginning, not the end.

Top Tags Applied by John
interpolation, arc length, curve, numbers, pchip
Files Tagged by John View all
Updated   File Tags Downloads
(last 30 days)
Comments Rating
06 Feb 2012 fminsearchbnd, fminsearchcon Bound constrained optimization using fminsearch Author: John D'Errico optimization, fminsearch, bounds, constraints, nonlinear, inequalities 292 46
  • 4.97368
5.0 | 38 ratings
01 Feb 2012 Screenshot interparc Distance based interpolation along a general curve in space Author: John D'Errico interpolation, interpolant, space curve, arc length, arc, distance 42 0
01 Feb 2012 Screenshot arclength Compute the arclength of a general curve in any number of dimensions Author: John D'Errico spline, arc length, curve, space curve, interpolation, arc 38 0
01 Feb 2012 Screenshot distance2curve Find the closest point on a (n-dimensional) curve to any given point or set of points Author: John D'Errico spline, pchip, interpolation, closest point, distance, curve 56 2
  • 5.0
5.0 | 1 rating
26 Jan 2012 Screenshot A suite of minimal bounding objects Suite of tools to compute minimal bounding circles, rectangles, triangles, spheres, incircles, etc. Author: John D'Errico minimal, minimum, area, perimeter, polygon, rectangle 63 0

Contact us at files@mathworks.com