@RM Wow, I tried to post that jar label peeling post on reddit too, but there was zero reaction to it. I guess you understand better how the system works!

So this thing I've been harrasing everybody about for the last couple days works. and it works really well.

nonlinear fitting like a baws

@ColinK IMO Mathematica's pretty good as far as nonlinear fitting is concerned. Or are you doing it in C with results you derived using Mathematica?

In Matlab with results derived in Mathematica

I'm using analytical expressions for the jacobian

@ColinK I see. Any special reason for using Mathematica and MATLAB together? MATLAB's symbolic capabilities not good enough yet? (I never used MATLAB at all myself, so I don't really know what it can do.)

so mma helped with the differentiation (these equations are gross) and helped me do some reformulating for numerical stability. The "normal" way of writing some of the jacobian terms have exp[-t] Erfc[-t] which explodes as t -> -Infinity

@OleksandrR I'm typically a Matlab user, and I don't do a lot of analytical math most of the time. I do a lot of statistics and numerical modeling. Matlab has a symbolic toolkit, but its not as good as Mathematica, not by a long shot

I could write this in C if I needed to, but it doesn't need to be that fast, and matlab is just easier. besides, I'm an expert Matlab programmer, and a mediocre C programmer.

@ColinK did you have any luck finding a CompexityFunction that adequately expresses numerical instability? I thought that was a very interesting question; not something I'd considered before.

I found something that does seem to at least encourage more stable solutions, but it's not perfect

it's also not very general. The ComplexityFunction needs to evaluate the expression at some typical value of the input parameters, so you couldn't use the same ComplexityFunction every time you want to look for numerically stable results

the CF explicitly substitutes and evaluates, like e /. {a -> 5, b-> 7} for example

grr, I cant make it format properly but you get the idea

Hm, okay. To make it more general, could you expand each expression as a series about any problematic values (poles, asymptotes, etc...) and check how fast it converged? The least complex would then be the expression with the most quickly converging series. (But, I haven't thought about how to implement that in practice, and it may be a bit tricky to get right.)

I might be wrong but I don't think that would work, for two reasons:

Well, that was not a good example was it? :) Apparently you can't mix fixed and proportional fonts in the same chat message.

@OleksandrR lol, I figured that was the problem

For example, Exp[-t] Erfc[-t] doesn't do anything special at t = -1000, but Exp[-t] gets HUGE and Efrc[-t] gets miniscule while their product is still small. That results in huge rounding error when you're using finite precision floating point

but there are no poles, zeros, asymptotes, etc. at that point

and, the way you want to reformulate it is Exp[-t + Log[Erfc[-t]]] which is exactly equivalent, but avoids the extreme intermediate values. But, since its exactly equivalent, all of its derivatives, and thus its series expansion, are identical at all points and to all orders

Very good points. I guess it's rather obvious that I hadn't thought very much about that. :)

It was a cool idea though. It might catch some cases, as there are, of course, numerical problems around poles, for example

the only real trouble is that they arent only at poles, but a lot of the certainly are

I think it should be possible to do this reliably and generally with Mathematica, but I think to do it right would take more than a one-line ComplexityFunction and Simplify[]. Probably a custom function that uses Level[] and looks for well defined situations that result in instability. My case (Huge * miniscule = small) is just one type of problematic situation in floating point math

Indeed. I think that is the main reason why Mathematica implements arbitrary precision arithmetic with adaptive precision control: if you get to a numerically unstable situation, just add more working precision. Probably easier than figuring out automatically what is or isn't numerically stable.

a good guide for thinking about it is that floating point problems crop up in a lot of the same situations that produce really bad signal to noise ratios in statistics. For example, another problem case is Huge - Huge = small. In the statistical situation, both huge values can be expected to have a variance on the order of their value, so the small value is much smaller but has a huge variance, ie. it's worthlessly noisy.

@OleksandrR yeah exactly, although you bring up a good point: Does Mathematica actually adjust it's precision in response to some detected instability condition? Or does it do things in some way that is simply immune to numerical issues all the time?

Hm, that's a good way of thinking of it actually. As an experimental physicist I'm quite accustomed to thinking about experimental errors but less so with reformulating expressions into numerically stable forms.

if the former is true, then Mathematica already has some built in aparatus for detecting instability

@ColinK it uses a very crude zeroth-order approximation of how much precision is carried through or lost in certain operations and adds more than enough working precision so that it won't present a problem (except in pathological cases).

@OleksandrR Same here. Undergrad B.S. Physics, Graduate M.S Optics. I'm not really trained for this sort of thing but I've picked it up for what I work on.

When it sees precision dropping on an intermediate result, it bumps up the working precision.

the statistics <-> numerical accuracy equivalency isn't perfect by the way, but it's a good guide

@OleksandrR Hmm, is that logic accessible in any way to the user?

In a sense yes; $MaxPrecision, $MinPrecision, $WorkingPrecision, Precision[.], and Accuracy[.] are all manifestations of it.

I wonder if I could use level with Precision[]et. al., and make a ComplexityFunction[] that minimizes the needed precision to express the equation

Actually, maybe I'm being a bit too uncharitable: the approximation is not zeroth-order in all cases; some functions know exactly how much precision they lose. But it does assume that precision on all results is uncorrelated.

still would need to provide typical values though

@OleksandrR Probably a fine assumption for almost all uses

it may not always be true but if it makes a difference you're in a really really bad situation anyway

I can't imagine how tough a situation you'd have to be in to be depending on the inaccuracy in one result to compensate for an opposing error in another result. Yikes.

@ColinK sure, although I'm so used to considering it in experimental error analysis that it feels odd to be able to make that assumption and not worry too much about it.

Numerically I mean. In actual science thats not uncommon, of course.

@ColinK maybe you can ask rcollyer about error cancellation in density functional theory... ;)

yeah, exactly, in actual real world science, thats very common. but when your "Result" is the result of a mathematical function in a computer, and not the measured result of some experiment...

I mean if your problem is so badly scaled that the results are all up against the precision limit of 64-bit IEEE floating point...?

Yeah, I understand. I certainly wouldn't want to be in that situation; in fact I think that's the main reason why libquadmath et al were invented.

yeah, sometimes it's unavoidable I suppose. Or at least convenient to not have to worry.

If oyu can get away with some algebraic manipulation that gives you a stable expression, its much better, and probably faster too. basic math functions have been so well optimized at this point, you may as well use them as much as possible

for example, this optimization, now that it doesn't have to do finite differences to get the jacobian, and the analytical jaacobian isn't corrupted by rounding error, is SO FAST

and absurdly robust too. I can't find any data that upsets it

@ColinK that's really good. I went the other way recently, and implemented a faster version (than what Mathematica provides) of a Nelder-Mead optimizer. But of course derivative-free methods are always going to need more function evaluations.

Yeah I was doing Nelder-Mead previously.

derivative methods let me make stronger statements about optimality and make it easier to compute confidence bounds, also it doesn't need as good a starting point

but honestly I use nelder-meade for most optimization problems, at least it's what I try first

and usually it works.

actually I've always been curious how N-M seems to be so incredibly good for so many situations. No other simple search/optimization algorithm I've seen is as robust and general. 90% of the time I just shove my problem into Matlabs N-M function and don't even adjust the default parameters of the solver, and out comes the answer

@ColinK and if not, you can try differential evolution, etc.

@ColinK I think you and a lot of others. There's no general convergence proof for N-M (and not through lack of trying), but in practice it works surprisingly well as long as you don't choose terrible starting values.

indeed. It's like magic

oh, Differential Evolution sounds a lot like a genetic algorithm, in concept at least

the "Global Optimizer" in some lens design software works that way

Oh, this is my only opportunity to go get lunch. Nice talking :)

@ColinK kind of. It's a stochastic minimizer but the biology analogy is much looser and parameters are real valued rather than binary coded. Actually it's more closely related to evolutionary strategies than GAs.

Okay, bye then! I should get back to work.

@OleksandrR That's a very good question, and I wish I knew what it was. In dft (density functional theory), there are three sources of inaccuracy: the limited basis size, the functional approximation to the exchange correlation energy (which combined with the "small" basis means convergence to a different ground state), and it is all cast as a non-linear minimization. Fun stuff, all around.

@rcollyer I guess in the usual case the error due to an incomplete basis set should be small (or can be made small) compared to the error in the functional. But that in itself contains a number of cancelling errors: I recall nobody thought that LDA would work at first, until it was tried, and actually it works fairly well.

@OleksandrR Truthfully, I don't look into the internals of the functionals if I can help it: that way lies madness. :) But, you are right the basis set size can be managed. For instance, I run several calculations at different basis sizes, and fit the energy to a decaying exponential. This gives me an error estimate for the energy. Truthfully, I don't know of anyone else who is as careful about that as I am.

Additionally, there's the whole minimization procedure itself, which is pretty nasty. But, something I would like to understand better.

@rcollyer I looked into it at one point (mainly because someone suggested that I use B3LYP for inorganic solids, which I strongly disagreed with). But yes, I agree it's complicated, and I forgot the relevant details now.

@OleksandrR I was only looking at the differences between LDA compositions (two different versions of Perdew's) and those were going beyond what I could digest comfortably without spending to much time on it.

@OleksandrR I guess I am missing something about DistributeDefinitions and compiled functions

what do you get when you evaluate ParallelEvaluate@HoldForm[testwm // Evaluate]?

@acl {testwm, testwm, ...}. I think your approach with OwnValues is perfectly valid too, but DistributeDefinitions has some weird behaviour. For example if you distribute the definitions before launching the kernels it's okay; if you do it afterwards, it's not.

I personally don't like the automatic distribution; I find it confusing and unpredictable. Obviously, when it behaves like this, it's even worse.

@OleksandrR yes I later noticed other strange things.

Well, I don't understand what it is doing most of the time. in "working" code, I try to avoid this kind of parallelization as much as I can.

in any case, looking forward to see some explanation of what's going on in an answer to the question

I also don't like the fact that they changed SharedDownValues to SetSharedFunction in 7. It's Mathematica-ly illiterate and confusing, IMO. As for this new strangeness in DistributeDefinitions, I guess I will have to have a look at it later.

@Szabolcs I posted that link before you on Reddit. It got me a bronze announcer badge (25 visits), not much, but better than your zero response ;-)

But you were quicker with volume rendering question, I see...

@szabolcs why did you delete your vertex indexing answer?

I'm looking at this question question, but I'm wondering how a connected graph can have sinks and/or sources.

@Heike it doesn't. The second requirement is rather redundant, I think

@SjoerdCdeVries You posted it in /mathematica/, but RM posted it in /programming/, and it got many more views there. I also tried to post in programming right after you did, but it had no effect

@Heike Connected = connected assuming undirected edges (even if the graph is directed). Then it can have sources and sinks. These two requirements are different from asking it to be strongly connected.

@SjoerdCdeVries I deleted because I realized he was asking about vertex labels (i.e. hwo they're labelled in a graphic), not names

@Szabolcs but don't the two together imply strong connectivity?

@RM try {1->2, 2->3, 3->1, 4->5, 5->6, 6->4, 3->4}

@Szabolcs right, and aren't the 2 possible subgraphs for his problem all strongly connected?

@RM You could select 1, 2, 6 ...

@Szabolcs you mean as an initial vertex list? shouldn't that result in an error?

ah, ok. I see what you mean. But since the edges are known to be directed, is it right to use the definition of connected for undirected edges?