16

Quasi-Newton methods construct an approximate Hessian for an arbitrary smooth objective function $f(x)$ using values of $\nabla f$ evaluated at the current and previous points. At each iteration of the method the quasi-Newton approximate Hessian is updated using the gradient evaluated at the latest iterate, $x^{(k)}$. These approximate Hessians aren't ...


11

Your example is a pretty good indication that the two derivatives (with respect to $x$ and with respect to $u$) do not commute :) (In fact, they're very different beasts -- one is a Fréchet derivative, the other a weak derivative.) Rather, you should consider $a(u) = (\partial_x u)^2$ as the composition $a = f\circ g$, of two functions $$ f: v(x)\mapsto v(...


10

According to Wolfram Alpha, $x^5+3(x-1)=0$ has no closed-form solution, so you can forget about a nice closed-form expression. :) I see nothing wrong with Newton's method; it should be quick and accurate, and with some analysis like the one you are sketching I think you can identify a safe starting point and prove global convergence. It might even be faster ...


8

Yes: See Higham's book "Accuracy And Stability of Numerical Algorithms", second edition, chapter 25: Nonlinear Systems and Newton's Method. In particular, see the section on the "limiting residual" in terms of a condition number for the Jacobian. It may well be that since your system is ill-conditioned, that you quickly hit the limiting residual and your ...


7

This is somehow unexpected, but my recent experience with solving a system of nonlinear equations is that treating them as the right hand side of a system of ordinary equations and then evolve the system with an ODE solver can be considerably faster than with the usual Newton-Raphson iteration. It sounds like you're doing some sort of pseudotransient ...


7

The issues you're running into now are not a failing of Newton-Raphson, but a question of coupling. You're doing iterated sequential coupling -- solving each equation sequentially and then iterating until (hopeful) convergence. No solver choice in place of NR is going to fix this lack of convergence, as long as you are doing iterated sequential coupling. ...


7

It's a bit easier to see if you write your equation in the a semi-discretised system of the form $u^{\prime}(t) = F(u(t))$ and with the application of the $\theta$-method and approximating $u^{\prime}(t) \approx (w^{n+1} - w^{n})/\tau$ this gives, $$w^{n+1} - w^{n} - (1-\theta) \tau F(w^n) - \theta\tau F(w^{n+1}) = 0$$ with unknown vector $w^{n+1}$ and ...


7

This is sometimes known as Buzano's inequality (http://www.jstor.org/stable/2159168). In general, if $\|x\|=1$, $P=xx^{\top}$ is a projection operator, so a simple application of Cauchy-Schwarz leads to $$ 2|b^{\top}xx^{\top}a|-|b^{\top}a| \leq |b^{\top}(2P-I)a| \leq \|a\|\,\|b\|, $$ which gives Buzano's inequality $$ |b^{\top}xx^{\top}a| \leq \frac{\|a\|\,\|...


7

A low rank update to an $n$ by $n$ matrix $A$ is an update of the form $A=A+UU^{T}$ where $U$ is matrix with $n$ rows but very few columns (typically just one or two.) If the matrix $UU^{T}$ has $k$ columns, then its rank is at most $k$, so this is a low-rank update to the $A$ matrix. Low rank updates are important in computational linear algebra and ...


7

Some of the confusion comes about because mathematicians often consider artificial functions that happen to have size $O(1)$ whereas practitioners use functions and variables that just happen to be whatever they are, often vastly different from one. One of the considerations that is important in practice is to make sure that the choice of physical units does ...


6

There are multiple ways of implementing Dirichlet boundary conditions when using the finite element method. For each unknown $i$ of the system belonging to the Dirichlet boundary, you can zero out row $i$ of the Jacobian and entry $i$ of the right-hand side. At that rate, you can reformulate the problem for just the unknowns in the interior, which reduces ...


6

Off the top of my head, there are a number of things that may be going wrong. You might want to verify the following Are the function values or derivatives you are computing numerically stable? Is the rounding error in their evaluation smooth? You can verify this visually by plotting an extremely small interval around a known root such that the function ...


6

A few specific examples: Quasi-Newton methods avoid computing second derivatives by using an approximation to the Hessian which is updated at each iteration by a low rank (rank one or rank two) update. This makes factoring the Hessian (or equivalently keeping the inverse in product form) easy to do computationally. There are also limited memory Quasi-...


6

The approximation error is $\mathcal{O}(\varepsilon)$, where $\varepsilon$ is the parameter in the finite difference calculation. It will be good to know how this approximation affects the number of CG iterations compared to having the complete Hessian? Are there theoretical/empirical understandings about this? Yes. In short, there are problems where CG ...


6

To calculate the reaction forces at a node, Abaqus (or any structural FE code) simply sums the internal forces for all elements attached to that node. The reaction forces are the negative of that sum. For an Abaqus user element, the internal forces for the element are returned from subroutine UEL in the RHS array. The returned stiffness matrix (Jacobian), ...


6

The reason is that GMRES can only be used for solving linear equations, i.e. equations of the form $Ax=b$, where $A$ is some matrix and $x,b$ are vectors. What GMRES does, essentially, is it approximates multiplication by the matrix $A^{-1}$ using a matrix polynomial of $A$. In this case (I assume) $f(y^{n+1},t)$ is not necessarily linear in the vector $y^{...


5

Your observed quadratic convergence indicates that the Jacobian is likely correct. Have you looked at the solutions for your under-resolved configurations? Galerkin optimality uses the operator norm, which contains only the symmetric part, thus the solution of the discrete system could be quite different from the projection of the exact solution. This ...


5

I'm not sure if this is the answer, but consider the case as $\alpha \rightarrow \infty$. In this case you're trying to solve $\iint_{\Omega} \cosh(P) dxdy = 0$ And since cosh looks like this: You can see that integrating it over any domain will never give you zero. If the diffusion of $P$ ($\nabla^2 P$) is not large enough to counteract this limit then ...


5

I suspect the main problem is the magnitude of the values. If you divide through by $\cosh$ to make all the numbers smaller, then ApproxFun doesn't seem to have a problem finding all the roots. function lhs(ε, α) A, B = sqrt((α-1)/(2ε)), sqrt((α+1)/(2ε)) return sqrt(α*α-1)*(cos(A) - 1/cosh(B)) - sin(A)*tanh(B) end u = Fun(x -> lhs(1e-3, x), [1.0,...


5

I don't know much about tracking implicit surfaces, so I'm just going to start with the optimization problem and go from there. The optimization problem is, at the core, nonlinear least squares, and can therefore be solved efficiently by the Gauss-Newton method, wherein at each step one linearizes the problem at the current guess, solves the linear least ...


5

The Netwon(-Raphson) method helps you finding the root(s) of a function $f(x)$, that is, those values of $x_0$ for which $f(x_0) = 0$. In your case, the unknown value is $c$, so we start by writing down the problem in terms of a function $f(c)$ that is zero when $c$ takes the "correct" value (which we call $c_0$). Since this looks like homework, I don't ...


5

Not sure where you get your equation for $\epsilon$, but ultimately your approximation for the Jacobian matvec operation is a finite difference approximation to the directed derivative of $F(\cdot)$. This means that you will want as accurate of an approximation to this directed derivative as possible while being sure to avoid numerical round-off issues. For ...


5

If you have a single vector equation $\vec{F}(\vec{x})=0$ then you solve it by representing that state vector $\vec{x}$ as a set of amplitudes $[x_0,x_1,...,x_{n-1}]$ after discretization by your favorite method (FD, FV, FEM, spectral); and we know how to solve it. If you also have a second equation $\vec{G}(\vec{y})=0$ then the full state vector is $[x_0,...


4

They're equivalent. Both imply some variation of the root-finding method by linearization.


4

Using a line search or trust region method along with Newton's method typically does the trick of finding the zero reliably. Check any advanced numerical methods or optimization book.


4

Can you bound your zero? If you can come up with a good heuristic for the approximate location of the zero and then compute bounds for it over an interval that is not too large, the performance of many methods is improved. Here's some Matlab code I use to perform basic root bracketing for monotonic functions. In terms of interval methods, @Pedro's ...


4

This remark is a bit too long for a comment, but regarding DanielShapero's remark with respect to changing the Jacobian matrix, for every degree of freedom $i$ corresponding to a Dirichlet boundary condition (BC), a common approach is to set $i$th row of the Jacobian matrix equal to the $i$th row of the identity matrix, meaning your (quasi-)Newton iterates ...


4

The most robust way of answering this is to benchmark it. Failing that, there are several things to note (roughly in order of importance). First, the most cheap floating-point operations on a modern CPU are addition and multiplication (both are equally fast; same as fused multiply-adds when available). Division is much slower (by a factor of ~20), ...


4

If you converge, you would expect the steps to get small. Ideally, a step $\delta x_k$ in an optimization algorithm would go from the current iterate $x_k$ to the exact solution $x^\ast$, so $\|\delta x_k\| \approx \|x_k-x^\ast\|$ which is, as you approach the solution, going to be small. Now you say that $\delta x_k$ is getting small even though the ...


4

It looks like this paper is combining Hessian-free with truncated Newton method. Yes, it is. ...the approach is referred to as Hessian-free method. That is because the Hessian is never computed explicitly. Instead, the product of the Hessian and a vector is obtained using finite difference approximation. Right. This step is analogous to the "Jacobian-...


Only top voted, non community-wiki answers of a minimum length are eligible