问:
Algorithm for integrating a 6D function in a Morse-Smale 3D cell
Lets say that one has a scalar field defined in 3D space for whose gradient he wants to find the Morse-Smale Complex for later performing an integration of several hexa-dimensional functions over attractor domains. The functions to be integrated are mostly concentrated around the attractor with relatively small contributions at the boundary. As well, the boundary can have a spiky shape. What approach (or algorithm or even software) would be appropriate to tackle this problem?

Currently the problem is divided in two steps: first the boundary of each attractor is found by testing in many points in space if following the gradient drives the search to other attractor or not. The initial points are chosen along rays starting at the attractor. Those rays point in the direction of the points generated by a Lebedev grid on a sphere of unit radius. So later the integration is separated in radial and angular variables. For each point in the radial quadrature (that goes no further than the point of the boundary found furthest from the attractor) the corresponding sphere is integrated checking for each angular point if it is inside or outside of the region of integration.

I am wondering if there is a more efficient approach around it that has enough resolution to treat the spiky boundaries and the maintains the benefits of the separation in angular and radial integration.


问:
Methods to compute specific eigenvector components for a tridiagonal matrix
I have an application that is somewhat similar to the situation of computing Gaussian quadrature nodes and weights: simply put, I need to compute the eigenvalues and the last two (normalized) components of the eigenvectors of a symmetric tridiagonal matrix.

I am aware of the Golub-Welsch scheme that uses the implicit QL/QR algorithm for Gaussian quadrature rule generation; the difference between it and my situation is that one uses the first component of the eigenvectors for Golub-Welsch.

I can modify the relevant LAPACK routines for my specific situation, but I was wondering if there were any new developments in this area. From my limited searching, it does not seem that MRRR or divide-and-conquer allow for computing only selected components.

Literature pointers would be extremely appreciated.

ContentLicense: CC BY-SA 4.0


问:
First order methods for a large scale semidefinite program
I am interested in solving the following semidefinite optimization problem:
\begin{equation}
\begin{split}
\underset{X,\lambda}{\rm maximize} \;\;\;\;&\lambda^Tc \\
&-\mathbb{I} \le X \le \mathbb{I} \\
&X - \sum_{k=1}^m \lambda_k A_k \ge 0,
\end{split}
\end{equation}
where $\lambda,c \in \mathbb{R}^m$ and $X,A_k$, for $k = 1,\ldots,m$, are symmetric $n\times n$ matrices, where $c$ and the $A_k$'s are fixed. Further, $\mathbb{I}$ denotes the $n \times n$ identity matrix. I would like to solve such optimizations for $n$ large ($\mathcal{O}(10^3)$ and larger) and $m \sim 100$.

Are there any algorithms that perform better than alternating projection for those dimensions? What is the state of the art algorithm for such a problem?

OwnerDisplayName: Marc


问:
Use of Morton Key to reduce number of grid points
I asked a question on Stack Overflow Performance Issue with VP Trees and Nearest Neighborsand I was not satisfied with the answer and so I thought I would reword my question for this site and post this question as a computational science question.

I am using the VPTree to store the coordinates of my sparse array as data is on a sphere/ellipsoid and I want to use great circle distances  or geodesics to do nearest neighbor searches as asked in a related question - Fast nearest neighbor search with latitude and longitude . The creation of the VPTree is itself not a problem as that happens very fast( I have 7000 points approximately). I create a regular grid in the following way and the minLat,maxLat,minLon,maxLon are obtained from the coarse data grid that contains the sparse matrix of observations.

```<code>latGrid = np.arange(minLat,maxLat,0.05)
lonGrid = np.arange(minLo,maxLo,0.05)
</code>```

So the total number of points in latGrid is 1800 and the total number of points in lonGrid is 7200. So the total number of grid points is 1800 * 7200 = 12960000

I am trying to do a Fixed radius nearest neighbor search and that is where I have my problem currently -

```<code>for grid_point in (grid_points):
    indices = tree.get_all_in_range(grid_point,r=4)
    args.append(indices)
</code>```

When I run this query against 12960000 grid points it takes 12 minutes and I have to run this 175 times.

I build the VPTree using the 7000 points of the coarse grid and I query it with 12960000 points of the regular grid.

I believe the number of grid points that is sent to the query

```<code>get_all_in_range
</code>```

can be reduced substantially as many of the grid points  are not in the range of the sparse matrix coordinates. The question is how  does one do an "intersection" of the grid points and the coarse matrix coordinates using a mathematical approach ? I have seen the answer to this question - Using Morton keys or Geohashes but I am not sure whether this would be applicable to my problem. At least I do not understand how to apply it to my problem. Similarly can a geodesic convex hull be used to determine the intersection ?


问:
Reweighted least squares factorization
This is a continuation of the question asked here. I want to solve numerous least squares systems of the form
$$
D_i A x \approx D_i b
$$
where $D_i$ are $m \times m$ diagonal matrices with positive elements on the diagonal, $A$ is $m \times n$ with $m > n$, $x$ is $n \times 1$, and $b$ is $m \times 1$.

The answer given in the previous question recommended to factor $A$ with an $LQ$ decomposition, noting that $D_i L$ would be lower triangular so that it wouldn't be necessary to factor $D_i A$ for each $D_i$.

However, it turns out that $L$ in an $LQ$ factorization is only lower triangular when $m \le n$ (under-determined). When $m > n$ (over-determined), $L$ is in fact lower trapezoidal (zero above the $m=n$ diagonal, see the LAPACK routine DGELQF). So $L$ would look like
$$
L = \pmatrix{ L_{11} \\ L_{21}}
$$
where $L_{11}$ is $n \times n$ lower triangular, and $L_{21}$ is $(m-n) \times n$ dense.

Therefore, when using an $LQ$ factorization, we would need to solve
$$
D_i L z \approx D_i b
$$
where $z = Q x$. Now multiplying a diagonal matrix by a trapezoidal gives a trapezoidal matrix, so the problem becomes solving the trapezoidal least squares problem
$$
(D_i L) z = y
$$
My question is how can this trapezoidal system be solved efficiently? It seems to me I'm back to square 1 because I would need to factor $D_i L$ with $QR$ or $SVD$ in order to solve this least squares system. I don't see any LAPACK/BLAS routines which can efficiently solve trapezoidal systems.

Title: Reweighted least squares factorization


问:
Harmonic average of Diffusion Tensors in Finite Volume Method
I want to implement a Bilinear Finite Volume discretisation of the anisotropic diffusion problem:

$$\frac{du}{dt} = \nabla \cdot (\textbf{D} \nabla u)$$

Both my degrees of freedom as well as the Diffusion Tensors are given on the vertices of my quadrilateral grid. I can, therefore, use bilinear test/trialfunctions on each cell to calculate the full gradient.

I am wondering and am deeply puzzled as to how I am to approximate the Diffusion tensor on an internal sub-face of my control volume.

In the 1D case (TPFA), the diffusion coefficient is modeled to be constant over one cell and the degree of freedom is thought to be the same everywhere within the cell. The flux over an interface is made consistent by using the harmonic average of the two given diffusivities:

$$f_{A->B}= \frac{2}{\frac{1}{D_A}+ \frac{1}{D_B}} ~\nabla_x u$$

Usually, the gradient is then approximated by the simple finite difference of the values modeled to be in the cell centers.
sketch to clarify:

So far so good. Now I am not in a setting where the Grid is aligned with the principal directions of the diffusion tensors in involved (not "K-orthogonal"). The Two-point flux approximation would, therefore, fail, because it can not capture the off-diagonal entries of the diffusion tensors involved, because only the x-derivative can be calculated.

In my 2D setting, I do reconstruct the full gradient. I want to calculate four sub-fluxes within my cell as indicated by the following sketch:

Here the black dots indicate where my degrees of freedom lie. The dashed line shows the area where I model the diffusion tensor to be constant (dual grid), the friendly green dot marks the location where I want to evaluate the bottom-most of the four sub-fluxes to properly assemble. Following the procedure from the 1D case, I could use the harmonic average of the diffusion tensors of the involved dual cells A,B. Let's call it

option 1:

$$D_{eff} = \big(\frac{\textbf{A}^{-1 } + \textbf{B}^{-1} }{2}\big)^{-1}$$

Here the "-1" in the exponents indicates the matrix inverse. What puzzles me so much is that by calculating the matrix inverse, the diagonal entries become coupled with the off-diagonal ones. To my understanding, these are different physical processes and I feel awkward to have, say, The component $D_{xx}$ depend on $D_{xy}$ etc.

The second option I see is to do a harmonic average component-wise, to have:

option 2

$$D_{eff} = \begin{bmatrix}
\frac{2}{\frac{1}{D^{A}_{xx}} + \frac{1}{D^{B}_{xx}}} &
\frac{2}{\frac{1}{D^{A}_{xy}} + \frac{1}{D^{B}_{xy}}} \\
\frac{2}{\frac{1}{D^{A}_{yx}} + \frac{1}{D^{B}_{yx}}} &
\frac{2}{\frac{1}{D^{A}_{yy}} + \frac{1}{D^{B}_{yy}}}
\end{bmatrix}$$

If i do it this way, the diagonal component $D_{xx}$ is only dependend on the corresponding entries in $D_A$ and $D_B$.

What is the correct way to do average the Diffusion Tensors at the sub-face marked with the green dot?

UPDATE:
Investigating option two I found that as the offdiagonal entries of the Diffusion Matrices may be negative, I face the situation where I take the harmonic average of a negative and a positive number which, according to wikipdeia, is not defined properly. Help appreciated.


问:
Inflow and outflow boundary conditions for advection-diffusion equation
I'm trying to solve this advection-diffusion equation (ADE):

$$\frac{\partial \phi}{\partial t} + \nabla \cdot (-D \nabla \phi + \mathbf{u} \phi) = 0$$

In fact, this ADE framework is coupled to a Navier-Stokes (NS) solver and I get velocity field ($\mathbf{u}$) from NS solver. My question is about ADE boundary conditions on inflow and outflow planes. I postulate that because the flow is directed from inflow to outflow obviously, I put zero concentration ($\phi = 0$) boundary condition on inlet plane because my source of species is located far from inflow boundary. Also, I assumed that the dominant transport mechanism at the outlet plane is convection and as a result of that normal diffusive flux will be vanished at the outlet plane ($-D \frac{\partial \phi}{\partial \vec{n}} = 0$).

I wanted to know how much these assumptions are true and could be justified physically? I should say that the outflow boundary condition is verified by comparing the simulation results with analytical solution of the Graetz problem. Any reference, suggestion, or idea is appreciated.


问:
The sign of Schrodinger equation
I have a question for the format of Schrodinger equation

$$\psi(x,t) = \int_0^\infty c_n e^{-iE_nt/\hbar} \psi_n(x)$$

Why do we have $i$ instead of $-i$?

Title: The sign of Schrodinger equation


问:
Maple: patmatch does not identify patterns inside diff() operator
I understand that the function applyrule uses patmatch to identify patterns and apply transformations accordingly.

Suppose, I want to write a rule for transforming spatial derivatives along $z-$direction of field components $E_i(x,y,z)$ to momentum space, so that:

$\partial_{z}(E_i(x,y,z)) \to -ikE_i(x,y,z)$.

I tried to define a rule as follows:

```<code>In: applyrule(diff(f::function,z)=-I*k*f,diff(E1(x, y, z), z))
Out: diff(E1(x, y, z), z)
</code>```

It does not work. The reason is, patmatch is unable to match the above pattern:

```<code>In: patmatch(diff(E1(x, y, z), z), diff(f::function, z))
Out: false
</code>```

On the contrary, if there is no diff() operator,

```<code>In: patmatch(E1(x, y, z), f::function)
Out: true
</code>```

As a workaround, I currently write the following command:

```<code>applyrule(diff(dE1(X), z) = -I*k*dE1(X),expression)
</code>```

for each field component in expression and it works fine.

But I would like to know an easier way to do this so that I don't have to write the rules for each component separately. Specifically, the inability of patmatch to identify pattern inside diff() is troubling. Any workaround or solution to resolve this would be helpful.


问:
Efficient numerical optimization of an "almost separable" function
I have come across an optimization problem with the following objective function:

$$f(x_0,y_0,z_0,x_1,y_1,z_1,...,x_N,y_N,z_N) = \sum_{i=0}^N f_i(x_i,y_i,z_i, \alpha(x_{i+1}-x_i) + \beta(y_{i+1}-y_i))$$

i.e. the objective function is the sum of the functions $f_i$ that only depend on three variables $x_i,y_i,z_i$ and on a linear combination of the difference to the direct neighbors (for $x_i$ and $y_i$). So far I have tired an NLCG algorithm, but convergence is very slow.

Is there a specialized solver that can exploit the structure of the optimization problem?

I have thought about introducing new variables $s_i = \alpha(x_{i+1}-x_i) + \beta(y_{i+1}-y_i)$ and writing
$$f(x_0,y_0,z_0,s_0,...,x_N,y_N,z_N,s_N) = \sum_{i=0}^N f_i(x_i,y_i,z_i,s_i)$$ because then the optimization could be conducted separately for each $f_i$ but I would have to introduce equality constraints to get back the original problem?

$f_i(x_i,y_i,z_i,s_i)$ is nonconvex (there is a $\cos(s_i)$ term in the function); however $f_i(x_i,y_i,z_i)$ is convex.


问:
Computing Small Eigenvalues with Sparse Symmetric Indefinite Mass Matrix
I want the eigenvalues of the following generalized eigenvalue problem:
$$ Av = \lambda M v $$
where

$A\in\mathbb{R}^{n\times n}$ is sparse, symmetric, and positive semi-definite

$M\in\mathbb{R}^{n\times n}$ is sparse, symmetric, and indefinite (i.e. neither positive nor negative definite)

I want the eigenvalues $\lambda$ with the smallest absolute value, i.e. $|\lambda|$ closest to zero.

My problem is that I am using scipy and in particular scipy.sparse.linalg.eigsh. (Note that I have to use Python, but not necessarily scipy, if someone knows a better library.) In this case, eigsh forces me to use the buckling option, since $M$ is indefinite. However, the buckling  mode sorts based on eigenvalues transformed via:
$$ \widetilde{\lambda} = \frac{\lambda}{\lambda - \sigma} $$
where $\sigma$ is a user-defined constant. Unfortunately, I can't figure out how to get the eigenvalues with smallest magnitude from this.
If I choose $\sigma$ to be large, then ask eigsh for the smallest magnitude (i.e. option SM), it is too slow. Slightly better is using the option BE, which returns values from both ends of the spectrum, with a small $\sigma$ close to zero. However, there's no reason to believe the smallest eigenvalues are equally distributed between positive and negative.

Does anyone have any advice for dealing with getting the smallest magnitude eigenvalues from a generalized symmetric eigenvalue problem where the mass matrix is indefinite?

(As a separate note, I have to add $\xi I$ to $A$, and then subtract $\xi$ from $\lambda$ at the end (for small $\xi$), so that it goes from positive semi-definite to positive definite, as required by buckling mode. I am aware that this is not an exact correction for the generalized eigenproblem as it is for the simpler case, but I could not find a better alternative so far, and it doesn't seem to be a big issue. I've asked a related question here, specific to this aspect.)

EDIT (091018): I needed some preliminary results and it seems like the following works, but certainly may not be optimal in terms of speed. Simply solve $M^{-1}Av=\lambda v$ instead using eigs. The trick to making it fast is using prefactorizations (i.e. sparse LUs) to make computing $M^{-1}A$ and $A^{-1}M$ fast. I'd love to know more about potentially better ways to do this!


问:
B-splines least squares with equality constraints
Can someone recommend the best way to solve a least squares fitting problem with B-splines, with additional equality constraints? I want to solve:
$$
\min_x || b - A x ||^2, \textrm{subject to: } C x = d
$$
where $A$ is $m \times n$ and $C$ is $p \times n$. Specifically I want to fit B-splines and force the first derivative of the resulting spline to take prescribed values at a certain number of points.

For the unconstrained problem, the standard approach seems to be to solve the normal equations
$$
A^T A x = A^T b
$$
because for a basis of B-splines of order $k$, the normal equations matrix  $A^T A$ is symmetric and banded, with lower bandwidth $k-1$. So a banded Cholesky decomposition can be applied to efficiently solve this system.

However I can't find much literature on dealing with constraints and B-spline fitting. For linear equality constrained systems, LAPACK uses the Generalized RQ decomposition, which does not require forming the normal equations matrix. I could of course do this as well, but the Generalized RQ algorithm won't then take advantage of the sparse structure of the matrix $A$. It might be possible to adapt the RQ algorithm to deal with the sparsity of $A$, but there don't appear to be any existing libraries which do this.

Just wondering if anyone has experience with this and/or recommendations?


问:
WENO methods: why the characteristic wise method resulting big errors?
I was doing my research/project using WENO as the limiter in finite volume methods to solve hyperbolic conservation law. I have no idea why the result in the characteristic wise method has a big error compared to component-wise. This is my current solution that I get from my solver, it shows that there is a numerical issue in characteristic wise reconstruction.

Velocity in $x$ direction: solution of sod shock at $t=0.1$:

My algorithm of characteristic wise reconstruction (I am using 5th order WENO). Let we evaluate cell $i$:

Using ROE averaging, get averaged value $ \bar q_{i+\frac{1}{2}} = f(q_{i},q_{i+1})$ and compute left and right eigenvector

Transform variable $q$ to characteristic variable $v$ using left eigenvector, because 5 points is needed so we need $ q_{j} = [q_{i-2},q_{i-1},q_{i},q_{i+1},q_{i+2}]$, and then $v_j = L(\bar q_{i+\frac{1}{2}})*q_{j}$

Perform WENO reconstruction

Transform back to physical space using $ q^{\pm}_{i\mp\frac{1}{2}} = R(\bar q_{i+\frac{1}{2}})*v^{\pm}_{i\mp\frac{1}{2}} $, so we get $\ q^{+}_{i-\frac{1}{2}}$ and $q^{-}_{i+\frac{1}{2}}$

While at cell interface interface $i+\frac{1}{2}:\ $ $q^{+}_{i+\frac{1}{2}}=q^{+}_{(i-\frac{1}{2})+ 1},\ q^{-}_{i+\frac{1}{2}}=q^{-}_{i+\frac{1}{2}}$.

Compute flux $ F_{i+\frac{1}{2}} = f(q^{-}_{i+\frac{1}{2}},q^{+}_{i+\frac{1}{2}})$

Is there any inconsistency in my algorithm?


问:
Difference between Dishonest Newton method and Very Dishonest Newton method
What is the difference between the Dishonest Newton method and the Very Dishonest Newton method? Is there a difference or do they mean the same thing?

I have tried searching for this on the internet but I can't find anything. I understand that in the ordinary ("honest") Newton method, the Jacobian is computed in each iteration.
If we reuse the same Jacobian for several iterations, this is referred to as the Dishonest Newton method (DN) by some sources and as the Very Dishonest Newton method (VDN) by others.

Intuitively, I would think that the Very Dishonest Newton method reuses the same Jacobian for more iterations than Dishonest Newton. Maybe in VDN the Jacobian is only computed once? But I have not been able to verify this.

I would prefer some credible sources with clear descriptions of both methods as I want to refer to them in my master's thesis.


问:
Neumann boundary conditions in the Maccormack scheme
I am trying to solve the viscous Burger equation
$$
\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} = \mu \frac{\partial^2 u}{\partial x^2}
$$
with Neumann boundary conditions. I am
using the Maccormack predictor-corrector scheme:
$$
u_j^\overline{{n+1}} = u_j^n - \frac{\Delta t}{\Delta x}(F_{j+1} - F_j^n) + r(u_{j+1}^n - 2u_{j}^n + u_{j-1}^n )\\
u_j^{n+1} = \frac{1}{2}\Big[ u_j^n + u_j^\overline{{n+1}}  - \frac{\Delta t}{\Delta x}(F_{j}^\overline{{n+1}} - F_{j-1}^\overline{{n+1}} ) + r(u_{j+1}^\overline{{n+1}}  - 2u_{j}^\overline{{n+1}}  + u_{j-1}^\overline{{n+1}} ) \Big]
$$
for $F=0.5u^2$. I understand that the diffusive term has been centre-differenced, and I believe I implement the boundary conditions in this term by solving for a ghost node. For example, if 0 is my first real node, I say that the ghost node at $j=-1$ should satisfy $u_{-1}=u_{1}$, and I impose
$$
u_0^\overline{{n+1}} = \,... \,+ r(2 u_{1}^n - 2u_{0}^n ).
$$
My questions:

Is this correct?

How do I implement the no-flux condition with the forwards (backwards) differenced terms in the predictor (corrector) step?

Simply setting everything involving $F$ to zero doesn't seem to work.


问:
Sensitivity Lagrangian solution general case
I have asked this question already on maths and mathoverflow.

Just a question about a literature reference. I am writing a paper for engineers.

Usually for the Lagrange multiplier problem

$$
\nabla f(x)+\lambda \nabla g(x)=0
$$

the sensitivity result that the multiplier $\lambda$ gives the sensitivity for changes in the constraint function is derived for the case $g(x)−h=0$ for varying $h$.
Is there somewhere a reference deriving this for $g(x;h)=0$?

I looked through the literature and asked one expert. No trace of such a proof. But there must be one, I am sure.

I can do it myself, however, I would like to avoid to reinvent the wheel.

I am only looking for a hint for a reference, no explanations how to do it.


问:
Combining fluid flow solver based on lattice Boltzmann method with a mechanical deformation solver based on finite element method
I'm thinking to couple my fluid flow solver based on lattice Boltzmann method with a mechanical deformation solver based on finite element method to take account for solid deformation in my models. In fact, I'm looking for a fluid-structure interaction framework, which its fluid solver was written based on lattice Boltzmann. I didn't find any opensource example by using this coupling idea. Can anyone suggest any opensource software or direction how it is possible to do it?</p>


问:
Existing software/scripts for spiral graphs?
I am looking for existing software or scripts to generate spiral graphs from cyclical (time) data, as presented by Webber and Muller.  The graphs shown in the paper look like a great means of visualizing and identifying periodicity in cyclical data, but I haven't been able to Google up any existing solutions for creating graphs like these.</p>

ContentLicense: CC BY-SA 3.0


问:
finding null space to a complex matrix
I need to solve the following equation:
$$
\begin{pmatrix}
\frac{\omega^2}{c^2}\varepsilon_x-\mu_z^{-1}k_y^2-\mu_y^{-1}k_z^2 & \mu_z^{-1}k_xk_y & \mu_y^{-1}k_xk_z\\
\mu_z^{-1}k_xk_y &\frac{\omega^2}{c^2}\varepsilon_y-\mu_z^{-1}k_x^2-\mu_x^{-1}k_z^2 &\mu_x^{-1}k_yk_z\\
\mu_y^{-1}k_xk_z & \mu_x^{-1}k_yk_z&\frac{\omega^2}{c^2}\varepsilon_z-\mu_y^{-1}k_x^2-\mu_x^{-1}k_y^2
\end{pmatrix}
\begin{pmatrix}E_x\\E_y\\E_z\end{pmatrix}=\begin{pmatrix}0\\0\\0\end{pmatrix}
$$

All parameters are given except $k_z$, $E_x$, $E_y$, $E_z$. I first find $k_z$ that makes the determinant of the matrix vanishes, then I find vector $E$ as the corresponding null space of the matrix. I have the following codes that generate the matrix and find the values of $k_z$:

```<code>import scipy as sp   # np.sqrt doesn't handle complex number?

def maxwell_matrix(k0,kx,ky,kz,eps,mu):
    kx2,ky2,kz2,k02 = kx**2,ky**2,kz**2,k0**2
    mat = sp.empty((3,3),dtype=complex)
    mat[0,0] = k02*eps[0] - ky2/mu[2] - kz2/mu[1]
    mat[1,1] = k02*eps[1] - kz2/mu[0] - kx2/mu[2]
    mat[2,2] = k02*eps[2] - kx2/mu[1] - ky2/mu[0]
    mat[0,1] = kx*ky/mu[2]
    mat[1,0] = mat[0,1]
    mat[0,2] = kx*kz/mu[1]
    mat[2,0] = mat[0,2]
    mat[1,2] = kz*ky/mu[0]
    mat[2,1] = mat[1,2]
    return mat

def compute_kz(k0,kx,ky,eps,mu):
    eps = sp.array(eps)
    mu = sp.array(mu)
    if len(eps.shape)==1:
        eps = eps[sp.newaxis,:]
    elif len(eps.shape)==0:
        eps = eps*sp.ones((1,3))
    if len(mu.shape)==1:
        mu = mu[sp.newaxis,:]
    elif len(mu.shape)==0:
        mu = mu*sp.ones((1,3))
    kx2, ky2, k02 = kx**2, ky**2, k0**2

    A = kx2*(eps[:,2]*mu[:,0]+eps[:,0]*mu[:,2])\
            +ky2*(eps[:,2]*mu[:,1]+eps[:,1]*mu[:,2])\
            -k02*eps[:,2]*mu[:,2]*(eps[:,1]*mu[:,0]+eps[:,0]*mu[:,1])
    B = (kx2*mu[:,0] + mu[:,1]*(ky2-k02*eps[:,2]*mu[:,0]))\
        *(kx2*eps[:,0] + eps[:,1]*(ky2-k02*eps[:,0]*mu[:,2]))

    epsz_muz = eps[:,2]*mu[:,2]
    C = sp.sqrt(A**2 - 4*epsz_muz*B)

kz1 = sp.sqrt(0.5*(-A+C)/epsz_muz)
kz2 = sp.sqrt(0.5*(-A-C)/epsz_muz)

return sp.transpose(sp.array([kz1,kz2,-kz1,-kz2]))
</code>```

The equations for calculating $k_z$ are found using Mathematica and it should be correct as I have tested it many times.

I have tried two methods for finding the null space. The first method is SVD decomposition, and the second one is to find the eigenvector with eigenvalue zero. The following code does this:

```<code>import numpy as np
import scipy.linalg as spl

k0 = 1.0
theta = 37*sp.pi/180  # just some random number
phi = 0*sp.pi/180
eps = (2*np.random.rand(3)+1) + 0.1j*np.random.rand(3)
mu = (np.random.rand(3)+1)

kx = k0*sp.sin(theta)*sp.cos(phi)
ky = k0*sp.sin(theta)*sp.sin(phi)

kz = compute_kz(k0,kx,ky,eps,mu)[0]

tol = 1e-9
for i in range(4):
    mat = maxwell_matrix(k0,kx,ky,kz[i],eps,mu)
    print('det:',spl.det(mat))   # this confirms kz is correct

    # SVD decomposition
    U, s, Vh = spl.svd(mat)
    idx = sp.argmin(sp.absolute(s))
    if sp.absolute(s[idx]) < tol:
        evect = Vh[idx]
        evect = evect/spl.norm(evect)
    else:
        raise ValueError("no solution found!")
    dot_prod = sp.dot(mat,evect)
    print('SVD vector:', evect)
    print('SVD dot prod: ({x[0]:.3e}  {x[1]:.3e}  {x[2]:.3e})    abs:{y:.3e}'.format(x=dot_prod, y=spl.norm(dot_prod)))

    # find eigenvector with zero eigenvalue
    w,v = sp.linalg.eig(mat)
    idx = sp.argmin(sp.absolute(w))
    if sp.absolute(w[idx]) < tol:
        evect = v[:,idx]
        evect = evect/spl.norm(evect)
    else:
        raise ValueError("no solution found!")
    dot_prod = sp.dot(mat,evect)
    print('method2 vector:', evect)
    print('method2 dot prod: ({x[0]:.3e}  {x[1]:.3e}  {x[2]:.3e})    abs:{y:.3e}'.format(x=dot_prod, y=spl.norm(dot_prod)))
    print()
</code>```

The strange thing is that when $\varepsilon$ is real, both methods seem to give the same answer, but when $\varepsilon$ is complex, then SVD decomposition seems to fail. It agrees with method 2 for one set of $k_z$, but for the other set of $k_z$ the computed dot product is clearly zero. What's wrong with SVD decomposition/my method for finding null space? What's the best way to find the null space of a $3\times 3$ matrix, which is in general complex?

Title: finding null space to a complex matrix


问:
Solving a coupled eigen value problem
I have the following problem:

$$\begin{bmatrix}A &B \\C& D\end{bmatrix}\begin{bmatrix}x\\y\end{bmatrix} = \begin{bmatrix}\lambda I_m & 0 \\ 0& \mu I_n\end{bmatrix}\begin{bmatrix}x \\y\end{bmatrix}
$$

where $A$, $B$, $C$, $D$ are matrices, $I_m,I_n$ are identity matrices of sizes $m,n$, $x,y$ are vectors and $\lambda$ and $\mu$ are scalars.

How does one solve such a problem?

Title: Solving a coupled eigen value problem


问:
How to numerically optimize affine transformations?
I need to optimize affine transformations for of a set of triangles using energy function based on the connectivity.

The energy of an edge $e_j$ between triangles $T_a, T_b$ is given by
$$
E_j = \left[(A_a v_{a,1} + t_a) - (A_b v_{1,b} + t_b) \right]+
\left[(A_a v_{a,2} + t_a) - (A_b v_{2,b} + t_b) \right]
$$

where $(A_a, t_a)$ and $(A_b, t_b)$ are affine transformations on the triangles adjacent to the edge and $v_1$ and $v_2$ are the vertices of the edge. The solution $v_i \equiv 0$ is prevented by an additional deformation energy on the triangles, like the conformal energy.

After the optimization, we want to have

$$v_{a,1}=v_{b,1}, \quad v_{a,2}=v_{b,2} \quad \Rightarrow E_j = 0$$

My problem is, that this is a non-linear problem, which is hard for solvers. I am using Google Ceres with line-search, but it often fails to find a valid line-search stepsize when optimizing

$$
\begin{pmatrix}
a_1 & a_2 \\
a_3 & a_4
\end{pmatrix}, \begin{pmatrix}t_1 \\ t_2\end{pmatrix}
$$

as a vector $(a_1, a_2, a_3, a_4, \dots, t_1, t_2, \dots)$.

For some transformations, there are more favorable parameterizations, e.g.

Rotation can be parameterizized as

$$
\begin{pmatrix}
\cos \alpha & -\sin\alpha \\
\sin \alpha & \cos \alpha
\end{pmatrix}, \begin{pmatrix}t_1 \\ t_2\end{pmatrix}
$$

optimizing the vector $(\alpha, \dots, t_1, t_2, \dots)$.

Shearing can be optimized using

$$
\begin{pmatrix}
1 & s_1 \\
s_2 & 1
\end{pmatrix}, \begin{pmatrix}t_1 \\ t_2\end{pmatrix}
$$

optimizing the vector $(s_1, s_2, \dots, t_1, t_2, \dots)$.

But for general deformations, I do not see which parameterization of the transformation matrix would fit and optimizing the coefficients directly probably lets the solver fail (or not converge) because it yields no good gradient.

So what is the best way to optimize affine transformations to find triangle deformations?

Here a graphic what should be optimized:

The final optimization objective is:

Given a set of Triangles

$$T=\{t_1=[v_1, v_2, v_3], t_2=[v_4, v_5, v_6], ßdots\}$$

without shared vertices in the set and the information about the shared vertices of common edges as a separate set

$$E=\{[a,b,c,d], \dots|a\sim c, b\sim d\}$$

where $\overline{ab}$ and $\overline{cd}$ are the same edges, when all triangles are re-assembled to a mesh.

Find affine transformations $A_i+t_i$ for all triangles, such that for two adjacent triangles $t_1=[a,b,x]$ and $t_2=[c,d,y]$ holds

$$
A_1 a + t_1 = A_2 c + t_2 \quad\wedge\quad A_1 c + t_1 = A_2 d + t_2
$$
and for all transformations $A_i$:
$$
\text{trace} \frac{(A_i^T+A_i)}{2 \det A_i} \to \min
$$

The first condition ensures that common edges are merged and the second condition minimizes the conformal deformation, what e.g. prevents the mesh form collapsing to a single vertex.

When I for example optimize the triangles of a cut cube, then the only transformations I need are rotation and translation to arrange the triangles in the plane. Optimzing $(\alpha, t_1, t_2, \dots)$ using line-search finds a good solution without problem.s

But for example a sphere requires quite a bit of deformation to be embedded in the plane. And to optimize the deformation needed to assemble the triangles into a planar patch, the optimization of the raw values $(a_1, a_2, a_3, a_4, t_1, t_2)$ seems not to work very well, especially line-search with Google ceres fails to find a valid stepsize for which the wolfe condition holds.

I suspect that in contrast to e.g. the gradient of the rotation, the gradient of this vector is not well-suited for finding the optimal transformation and it would need another local parameterization, which better suits the search space, just like the rotation/translation for rigid embeddings.


问:
Imbalance of variables in Mixing Newton's method and Linear solver for a Non-linear system
Problem

Solving a non-linear system of equations.

The number of variables is the same as the number of equations.

When I fix a set of variables (say $\vec{y}$) and keep another set free (say $\vec{x}$), the system becomes an under-determined, dense, and linear system of the subset of variables
$$A(\vec{y})\vec{x} = \vec{b}(\vec{y})\tag{1}\label{system1}$$

where $A(\vec{y})$ is a dense matrix, and $\vec{b}(\vec{y})$ is a dense vector). Let's call this sub-system $\eqref{system1}$ as system 1.

When I fix $\vec{x}$ and keep $\vec{y}$ free, the system becomes an over-determined, sparse, and non-linear system of the subset of variables
$$F(\vec{x}, \vec{y}) = 0\tag{2}\label{system2}$$

The Jacobian $J(\vec{x}, \vec{y})$ is available in a closed form. Let's call this sub-system $\eqref{system2}$ as system 2.

About half of the equations in $\eqref{system2}$ are equality constraints that are linear in terms of $\vec{y}$. Each of the constraints is quite sparse and involves only about 5% of all variables.

Can I solve it with the following?

Algorithm 1

Initialize $\vec{x} = \vec{x}_0$, and $\vec{y} = \vec{y}_0$

Fix $\vec{y}_{n - 1}$, solve $A(\vec{y}_{n-1})\vec{x}_{n} = \vec{b}(\vec{y}_{n-1})$ for one of all the possible $\vec{x}_{n}$ because this system is under-determined.

Fix $\vec{x}_{n - 1}$. Perform one iteration of Newton's method for solving $F(\vec{x}_{n-1}, \vec{y}_{n}) = 0$ for $\vec{y}_{n}$.

If not converged, go to step 2.

Algorithm 2

If I replace step 2 in Algorithm 1 by a Newton's method iteration for solving system 1 $\eqref{system1}$, then I guess the steps become a block Newton's method.

Question

But I don't know if these two algorithms can work because system 1 $\eqref{system1}$ is under-determined and system 2 $\eqref{system2}$ is over-determined.

Can this work?

Related

A. Yeckel, L. Lun, and J. J. Derby, "An approximate block Newton method for coupled iterations of nonlinear solvers: Theory and conjugate heat transfer applications," J. Comp. Phys., vol. 228, no. 23, pp. 8566–8588, Dec. 2009

B. Polyak and A. Tremba, "Solving underdetermined nonlinear equations
by Newton-like method"
Solving over-determined non-linear equations with Newton's method and Moore-Penrose inverse is fine.

V. Meñkov. "Solving block linear system with low-rank off diagonal blocks is easily parallelizable"
Block triangular form allows some blocks to be over-determined while other blocks to be under-determined. Newton's method is equivalent to solving the system linearized at current approximation of solution. The question is equivalent to if it is possible to solve a system of linear equations by alternatively performing iterations of Newton's method on two groups of variables.

Śmietański, M. J. "A generalized Jacobian based Newton method for semismooth block-triangular system of equations" J. Comput. Appl. Math.. 2007.
Solving block triagular system by newton iteration on each block can work. "We prove locally superlinear convergence of parameterized version method that is well defined even when the generalized Jacobian is singular."


问:
Nonlinear functional optimization in radial coordinates
I am currently implementing classical density functional theory for a radially symmetric system.
In mathematical terms, I am searching for a function $f(r)$ that minimizes a functional $\Omega[f]$.
The simplest part of $\Omega$ looks like this: $\int\limits_0^R r^2 f(r)(\ln(f(r))-1)\mathrm{d}r$, but there are others that involve convolutions of $f$.
I am trying to minimize a numerical implementation of $\Omega$ using nlopt, which works fine in planar coordinates for the same functional.

$f(r)$ is represented as a vector $f_i$ which represents samples of $f$ at evenly spaced points $r_i$. The points start at $Δr$ and go up to $R$. $f(0)$ is irrelevant. There is a constant boundary condition $f(r>R)=f_∞$.

However, in radial coordinates, the convergence for small $r$ is way worse than for large $r$, because of the factor $r^2$ (the Jacobian determinant). If I only cared about the value of $\Omega$, this would not be a problem, but in this case, $\Omega$ is only a vehicle to obtain $f$.

Are there libraries that are better suited for this kind of problem? Or is there a different way to fix the problem and keep using nlopt?


问:
How to solve an implicit ODE with forward Euler?
Consider the implicit ODE

$$
M(y)\dot{y} = F(t,y)
$$

If $M$ is non-singular for all $y$

How to use the forward-Euler method to numerically solve for $y$ without inverting $M(y)$?

I only came out following method:

Forward Euler, $y_{n+1} = y_n + h  \dot{y_n}$

Multiply by $M(y_n)$ to obtain $M(y_n) y_{n+1} = M(y_n) y_n + h M(y_n) \dot{y_n}$

$M(y_n) y_{n+1} = M(y_n) y_n + h F(t, y_n) = \mathrm{RHS}$

Solve $M(y_n) y_{n+1} = \mathrm{RHS}$ by iterative matrix solver from start point $y_n$ and get $y_{n+1}$.

But I think the procedure is not efficient since it needs to iterative solve $M(y_n) y_{n+1} = \mathrm{RHS}$.

Could anyone tell me a practical or more efficient method (need to use forward-Euler)?


问:
Logging vs outputs in iterative optimisation
I'm coding an iterative algorithm of constrained continuous optimisation. An augmented Lagrangian algorithm (outer) calls a bound-constrained L-BFGS-B algorithm (inner), which calls a line search algorithm. Optimisation problems have between 1 and 10,000 bound-constrained variables and 100–200 nonlinear constraints.

The program outputs a solution vector (if found),  Lagrangian multipliers and termination reasons. Also, the output is a list of iterates, objective values and possibly Lagrangian multipliers per each inner and outer iteration.

For problematic cases, I'd like to provide diagnostics info: warnings, numerical issues (matrix ill-conditioning, NaNs, etc.) and so on.

I'm using Scala to adapt some algorithms from the Breeze library, and also add a Lagrangian algorithm. Some algorithms in that library have logging statements everywhere, whereas others have very little logging.

Questions:

What information is normally logged and what is output in numerical optimisation algorithms?

How do I decide between logging vs outputting? Obviously, there's a trade-off between outputs complexity and the "intrusiveness" of logging.

Will it make sense creating lean algorithm classes with diagnostics sub-classes containing logging?

Other advice, e.g., design tips, references and open-source code examples, are very welcome too.


问:
Efects from the boundary in advection equation
I am implementing the advection equation $u_x+(1/c)u_t=0$ following a Crank-Nicholson finite difference scheme. The equation for this is
\begin{eqnarray*}
  -\frac{\gamma}{4} w_{n-3 j+1} + w_{n-2 j+1} + \frac{\gamma}{4} w_{n-1 j+1}
  = \frac{\gamma}{4} w_{n-3 j} + w_{n-2 j} - \frac{\gamma}{4} w_{n-1 j}.
\end{eqnarray*}

In matrix form it can be written as:

\begin{eqnarray*}
  A \bf{w}_{j+1} = \it B {\bf{w}_j } + \bf{b}
\end{eqnarray*}
where

\begin{eqnarray*}
A =
\begin{pmatrix}
  1 & \frac{\gamma}{4} & 0 & 0 & \cdots & 0 \\
  -\frac{\gamma}{4} & 1 & \frac{\gamma}{4} & 0 & \cdots & 0 \\
  0 & \ddots & \ddots & \ddots & 0 & 0  \\
  \vdots & \vdots & \ddots & \ddots & \ddots & \vdots  \\
  0 & 0 & \cdots &  -\frac{\gamma}{4} & 1 & \frac{\gamma}{4} \\
  0 & 0 & \cdots & \cdots & -\frac{\gamma}{4} & 1
\end{pmatrix}
\; , \;
\bf{w}_{j+1} =
\begin{pmatrix}
  w_{1 j+1} \\
  w_{2 j+1} \\
  \vdots \\
  \vdots \\
  \vdots \\
  w_{n-2 j+1} \\
\end{pmatrix}
\end{eqnarray*}

\begin{eqnarray*}
B =
\begin{pmatrix}
  1 & -\frac{\gamma}{4} & 0 & 0 & \cdots & 0 \\
  \frac{\gamma}{4} & 1 & -\frac{\gamma}{4} & 0 & \cdots & 0 \\
  0 & \ddots & \ddots & \ddots & 0 & 0  \\
  \vdots & \vdots & \ddots & \ddots & \ddots & \vdots  \\
  0 & 0 & \cdots &  \frac{\gamma}{4} & 1 & -\frac{\gamma}{4} \\
  0 & 0 & \cdots & \cdots & \frac{\gamma}{4} & 1
\end{pmatrix}
\; , \;
\bf{w}_{j} =
\begin{pmatrix}
  w_{1 j} \\
  w_{2 j} \\
  \vdots \\
  \vdots \\
  \vdots \\
  w_{n-2 j} \\
\end{pmatrix}
\; , \;
\bf{b} =
\begin{pmatrix}
  \frac{\gamma}{4} w_{0 j+1} + \frac{\gamma}{4} w_{0j}  \\
   0 \\
   \vdots \\
   0 \\
   -\frac{\gamma}{4} w_{n-1 j+1}  - \frac{\gamma}{4} w_{n-1 j}
\end{pmatrix}
\end{eqnarray*}

As an initial condition I have a Gaussian pulse $u(x,0)=\mathrm{e}^{-400(x-1/2)^2}$. The spacial domain is $x \in [0,1]$. Constant velocity $c=1$, $dt=0.001$, $dx=0.002$ and $\gamma=c dt/dx$.

My implementation is in Python and I can provide it upon request.

I attach a few screen captures of the "wavefield".

Can someone please help me understand what is going on after the pulse hits the left boundary? I am using Dirichlet boundary conditions with 0 on the boundary.

Thanks.

*****Update*****

I found someone else asking the same question and some interesting answers which I am still trying to digest.  Still I would like to have more discussion on this issue.

Here is the link on this


问:
Order of a principal term
In Yurii Nesterov's Introductory Lectures on Convex Optimization, there is a bound for the total number of iterations for some process. See page 109:

$$\left[\frac{1}{\ln(2(1-\kappa))} \ln\frac{t_0-t^*}{(1-\kappa) \epsilon}+2\right]\cdot \left[1+\sqrt\frac{L}{\mu}\ln\frac{2(L-\mu)}{\kappa \mu}\right]\\+\sqrt\frac{L}{\mu}\cdot\ln\left(\frac{1}{\epsilon}\max_{1\leq i \leq m}\{f_0(x_0)-t_0; f_i(x_0)\}\right)\label{eq1}\tag{1}$$

Then, the principal term in the above estimate is of the order
$$\ln\frac{t_0-t^*}{\epsilon} \sqrt\frac{L}{\mu} \ln\frac{L}{\mu} \label{eq2}\tag{2}$$

How did we arrive at statement $\eqref{eq2}$? Is that true the second term $\sqrt\frac{L}{\mu}\cdot\ln\big(\frac{1}{\epsilon}\max_{1\leq i \leq m}\{f_0(x_0)-t_0; f_i(x_0)\}\big)$ in $\eqref{eq1}$ is eliminated? I would appreciate any advice here.

Title: Order of a principal term


问:
$L^2$ norm error estimates of conforming FEM about Poisson’s equation with mixed boundary conditions
Consider Poisson’s equation

$$- \Delta u = f{\rm\qquad{     in  }}\;\Omega $$

with following mixed boundary cconditions

$$u = g{\rm\qquad{ on }}\;\Gamma  \subset \partial \Omega $$
$$\frac{{\partial u}}{{\partial n}} = h{\qquad\rm{        on  }}\;\partial \Omega \backslash \Gamma$$

where $\frac{{\partial u}}{{\partial n}}$ denotes the derivative of $u$ in the direction normal to the boundary, $\partial \Omega$. When $g=0$ and $h=0$, the problem becomes homogeneous boundary condition.

I already know the Aubin-Nitsche trick to get error estimates for $u - {u_h}$ in the $L^2$ norm with homogeneous boundary condition. I learned the duality argument from Brenner and Scott's book The mathematical theory of finite element methods. But when the boundary condition is inhomogeneous, it is difficult for me to get error estimates for $u - {u_h}$ in the $L^2$ norm using duality argument. Although the book says,

...inhomogeneous boundary conditions are easily treated' and the proof is analogous

I encountered a problem arising from the Neumann boundary condition. When the boundary condition is homogeneous, I can get the dual problem

$$a\left( {v,w} \right) = \left( {u - {u_h},v} \right)\qquad\forall v \in H_0^1\left( \Omega  \right)$$

and the elliptic regularity estimates

$${\left| w \right|_{{H^2}}} \le C{\left\| {u - {u_h}} \right\|_{{L^2}}}$$

which is crucial to the proof. In the case of inhomogeneous boundary condition, however, I get the dual problem

$$a\left( {v,w} \right) = \left( {r - {r_h},v} \right) + {\left( {h,v} \right)_{\partial \Omega \backslash \Gamma }}\qquad\forall v \in H_0^1\left( \Omega  \right)$$

and the elliptic regularity estimates

$${\left| w \right|_{{H^2}}} \le C\left( {{{\left\| {r - {r_h}} \right\|}_{{L^2}}} + {{\left\| h \right\|}_{{H^{\frac{1}{2}}}}}} \right)$$

which $r = u - \tilde u$, $\tilde u \in {H^1}\left( \Omega\right)$ and ${\left.{\tilde u} \right|_\Gamma } = g$. In the following steps of proof, I try to estimate ${\left( {h,v} \right)_{\partial \Omega \backslash \Gamma }}$ and ${{{\left\| h \right\|}_{{H^{\frac{1}{2}}}}}}$, but I didn't get the result.

Questions:

Is my idea about the proof of the inhomogeneous boundary condition, such as the dual problem and the elliptic regularity estimate which I get, correct?

If my idea about the proof is correct, how can I continue to prove the error estimates for $u - {u_h}$ in the $L^2$ norm(my idea to estimate ${\left( {h,v} \right)_{\partial \Omega \backslash \Gamma }}$ and ${{{\left\| h \right\|}_{{H^{\frac{1}{2}}}}}}$)?

If my idea about the proof is wrong, how to obtain the error estimates for $u - {u_h}$ in the $L^2$ norm(the case of inhomogeneous boundary condition)?Please tell me the proof or the main idea of the proof.


问:
Randomized Submatrix of a Sparse Matrix
I have a sparse square matrix $A$ with size $n \times n$ and number of nonzero entries $nnz$.
The goal is making a sub-matrix $B$ with $s$ nonzeros which are randomly chosen from $A$. Duplicates are ok, so we need at most $s$ nonzeros in $B$.

The probability of entry $A_{kl}$ being chosen is $p_{ij} = \frac{A_{ij}^2}{\sum_{k, l} A_{kl}^2}$, in which, $k$ and $l$ go through all the matrix entries. I have this algorithm in mind so far:

generate a random number $rand$ such that $0\le rand \le 1$.

if $(rand < p_{11}) \implies$ add $A_{11}$ to $B$.

go to the next entry of $A$ and repeat from step $1$ until $s$ entries are chosen.

If I go through all entries of $A$ once, on average how many entries will be chosen? In other words, how many times I should go through entries of $A$ to have $s$ entries chosen?

I want to change the probability $p_{ij}$ so by going through entries of $A$ once, all $s$ entries are chosen.

Title: Randomized Submatrix of a Sparse Matrix


问:
Detecting structures (blood vessels) in MR images
What I have are 3D grayscale images of brains where dark pixels indicate blood vessels (see below, red arrows mark blood vessels). The objective is to obtain some descriptive statistics about them, such as number of vessels in the brain or volume fraction of the brain occupied by vessels.  Currently, some poor student would either manually count all the vessels in the brain, or we would just look at the histogram over a certain volume, to determine the volume fraction of blood vessels.

This is a bit a general question, but what approach would you take to automatically detect these structures? Scikit image seems to have some promising methods for edge detection. The algorithm should work on a 3D matrix since one vessel can span multiple image slices. The grayscale values could also be used to improve the resolution: A voxel that is completely filled by a blood vessel has the lowest signal value, while a voxel that is only partially occupied by a vessel has a slightly higher signal value.


问:
Solving electrostatics Poisson equation with Intel MKL routines
I am trying to solve  the 3D Poisson equation

$$\nabla\cdot(\epsilon(\mathbf{r})\nabla) u(\mathbf{r}) = f(\mathbf{r})$$

I notice intel advertises routines that appear to solve

$$\nabla^2u(\mathbf{r}) = f(\mathbf{r})$$

Does anyone know if I can find similar such routines for solving the more general electrostatics Poisson equation above?


问:
Where can I find sample data for large linear programming optimization problems?
I am doing a comparison of different algebraic modeling langues (AMPL, AIMMS, GAMS, Pyomo) in both theoretical and practical terms. As a practical experiment I am trying to measure problem model instance creation time. I.e. taking optimization problem defined in one of the above languages I am trying to load the model and it export/convert it to LP solver input format. Model instance creation time is being measured and model instance characteristics are being examined.

So far I have only used sample problems provided by the authors of AMLs which are quite small. In order to have a meaningful benchmark I would like to find test data to build much larger problem model instance.

Maybe someone could guide me where I could find larger data sets for the models of linear programming problems  defined in the AMLs mentioned above?


问:
Roe Riemann solver for perfect gas mixture
I have working program for solving one-component 1D Euler equations with Roe's approximate Riemann solver constructed according to this pdf. My implementation of the algorithm is as follows ($\rho$ is density, $u$ is velocity, $h=\frac{e+p}{\rho}$ is total specific enthalpy, $c$ is speed of sound):

First, we compute Roe average values between the states:

$$
\rho = \sqrt{\rho_l\rho_r} \\
u = \frac{\sqrt{\rho_l}u_l + \sqrt{\rho_r}u_r}{\sqrt{\rho_l}+\sqrt{\rho_r}} \\
h = \frac{\sqrt{\rho_l}h_l + \sqrt{\rho_r}h_r}{\sqrt{\rho_l}+\sqrt{\rho_r}} \\
c = \sqrt{(\gamma-1)(h-u^2/2)}
$$

Then the right eigenvectors of the Jacobian matrix:

$$
\left(
\begin{matrix}
r_0 & r_1 & r_2 \\
\end{matrix}
\right)
=
\left(
\begin{matrix}
1 & 1 & 1 \\
u-c & u & u+c \\
h-uc & u^2/2 & h+uc
\end{matrix}
\right)
$$

Absolute values of wave speeds:
$$
|\lambda_0| = |u-c|, \; |\lambda_1| = |u|, \; |\lambda_2| = |u+c|
$$
and apply some entropy fix to them (not relevant here).

Wave amplitudes (coefficients of the jump between states decomposition by eigenvectors):
$$
\alpha_0 = \frac{(p_r-p_l) - \rho c (u_r-u_l)}{2c^2} \\
\alpha_1 = (\rho_r - \rho_l) - \frac{(p_r-p_l)}{c^2} \\
\alpha_2 = \frac{(p_r-p_l) + \rho c (u_r-u_l)}{2c^2} \\
$$

Flux values of the states:
$$
F_l = \left(\begin{array}{c} \rho_l u_l \\ p_l+\rho_l u_l^2 \\ (e_l+p_l)u_l \end{array}\right), \quad F_r = \left(\begin{array}{c} \rho_r u_r \\ p_r+\rho_r u_r^2 \\ (e_r+p_r)u_r \end{array}\right)
$$

And finally flux values of the middle state:
$$
F = \frac{1}{2}(F_l+F_r) - \frac{1}{2} \sum\limits_{p=0}^2 |\lambda_p| \alpha_p r_p
$$

The problem is how to generalize this algorithm for equilibrium multi-component perfect gas mixture?

I'm already able to calculate the eigenvectors for this case, here is an example for 3 components, corresponding to $(u,u,u,u+a,u-a)$ eigenvalues:
$$
\left(
\begin{matrix}
r_0 & r_1 & r_2 & r_3 & r_4 \\
\end{matrix}
\right)
=
\left(
\begin{matrix}
1 & 0 & 0 & \frac{\rho_0}{\rho} & \frac{\rho_0}{\rho} \\
0 & 1 & 0 & \frac{\rho_1}{\rho} & \frac{\rho_1}{\rho} \\
0 & 0 & 1 & \frac{\rho_2}{\rho} & \frac{\rho_2}{\rho} \\
u & u & u & u+c & u-c \\
u^2-\frac{1}{\gamma-1}\frac{\partial p}{\partial \rho_0} & u^2-\frac{1}{\gamma-1}\frac{\partial p}{\partial \rho_1} & u^2-\frac{1}{\gamma-1}\frac{\partial p}{\partial \rho_2} & h+uc & h-uc
\end{matrix}
\right)
$$
where
$$
\frac{\partial p}{\partial \rho_i} = \frac{R_0 T}{\mu_i} - (\gamma-1)(h_{0i}+c_{vi}T - u^2/2)
$$
and $T$ is temperature, $R_0$ is universal gas constant, $\mu_i, h_{0i}, c_{vi}$ are molar mass, heat of formation and heat capacity of $i$-th gas component.

These formulas are derived from this paper by discarding of vibrational terms.

The questions are:

How to correctly apply Roe averaging to gas mixture? I.e. how to calculate $\rho_i, u, \gamma, T, h, c$ for the matrix above?

How to calculate $\alpha_i$ for this case?


问:
why I cannot find explicit finite difference for elliptic equation
Let us think on the Poisson equation $\nabla^2 u(\bf{x})=\rho(x)$ with Neumann boundary conditions, with $\bf{x}=\it (x,y)$ in 2D.

Here is a stencil with central differences in both $x$ and $y$ (forgive my label $t$ on the vertical, it is supposed to be $y$)  directions.

with $\beta=\Delta x/\Delta y$. More specifically we have the
equation:

\begin{eqnarray*}
  w_{i+1 j} -2 w_{ij} \left (  1 + \beta^2 \right)
  + w_{i-1 j} + \beta^2 (w_{i j+1} + w_{i j-1})
  = \Delta x^2 \rho(x,y)
\end{eqnarray*}

We can get a similar stencil in the (hyerbolic) wave equation in 1D where
time is the $y$ (vertical axis). We can always find, recursively
$w_{i j+1}$ in terms of all other $w_{lm}$ where $l=i,i+1,i-1$ and
$m=j,j-1$. The points 1 cell away from the boundary can be computed with
a ghost boundary condition (assuming Neumann BC with derivative equal to 0).

I have been searching for hours and all literature shows me an implicit method where I need to invert a huge $(n-1 \times m-1)$ matrix. Is there a specific reason why the elliptic equation is not solved as we do the 1D wave explicit solution centered in space and time?

Or...if it could be solved, is there a place to look for that solution? Thanks.


问:
Good C++ optimization library for BFGS
To implement maximum likelihood estimators, I am looking for a good C++ optimization library that plays nicely with Eigen's matrix objects. Eigen has some capability of interfacing of its own but if anyone here has experience of using Eigen with an optimizer library in tandem, that would be great!

I am aware of similar questions being asked here or here but those are specific for derivative-free methods.

I think a 'good' library would be one that plays nicely with other libraries such as Eigen or Armadillo. So far, I have found dlibnet and nlopt, which seem promising.

Ideally, the answer to this question could be a longer one with having the pros and cons of each library listed as well. Thank you so much!

Update: I have decided to go with CppNumericalSolvers for now. There was also Ceres (linking to the relevant bit) but what put me off that was the emphasis on non-linear least squares, which requires the objective function to be of a certain format, which for more complex statistical models cannot always be guaranteed. Even iterated generalised least squares which can achieve certain maximum likelihood estimators would fail this format.

Moreover, CppNumericalSolvers supports bindings to TensorFlow, Eigen, and MATLAB. Note that both Ceres and CppNumericalSolvers ask you to formulate your function as a class rather than a lambda function, which may put off some beginners.

Finally, there is ROOT which is an environment for C++ with tons of routines where C++ source files can be loaded as 'macros.' Elementary BFGS optimizers exist with plenty of examples such as here. However, this is an interpreted environment.

I am still hoping someone might answer my questions with some actual experience using a library.

Title: Good C++ optimization library for BFGS


问:
Uniformly sample a point per polytope

I want to uniformly sample a point within each of $10^5$ convex polytopes in each iteration of a solver.

The polytopes in one iteration are completely different from the polytopes in another iteration.

The dimension of the space is about 15.

The vertices of each polytopes are known.

I only want one point per polytope.

Hit and run sampling would take quite a lot of steps to converge to a uniform distribution.

Triangulating each polytope to simplex can work. But doing that in high dimension seems tedious (at a quick glance.)

Are there other methods to do this?

Thanks.

Title: Uniformly sample a point per polytope


问:
Solving Compressible Euler in Primitive Variables with Galerkin P1 FEM
I have implemented a small compressible Euler solver, discretizing in primitive variables (rho, u, v, p) with standard Galerkin FEM P1 triangular elements, and mixed isotropic and anisotropic/streamline diffusion for stabilization.

I have verified the solver on several test cases with shock tube, double shocks, and flow over bumps, so it seems to work fine except for some oscillations (as the artificial diffusion isn't shock capturing or monoticity preserving).

Now I'm trying to validate a simple 2D static oblique shock test case for a unit square which should be simple and has an analytical solution. With top/left boundary conditions Ma_in = 2, rho_in = 1, u_in = cos(10 deg), v_in = sin(10 deg), p_in = 0.1786, bottom slip (v = 0), the analytic solution should be a shock at 29.3 degrees and Ma = 1.64 at the outlet (right boundary). My solver converges just fine and resolves the shock, however the solution is inverted after the shock, Ma 2.4, density 0.65 instead of 1.46 etc.

I cannot seem to find the error, especially as the dual shock problem solves fine. Could this be to using non-conservative/primitive variables? Any ideas at all would be helpful.


问:
Thermal stress to displacement through finite volume method
I am attempting to solve for a displacement field where I know the thermal stresses on my discretized domain, which consists of hex cells.

My first question is: is easier way to solve for the displacements besides integrating the linear elasticity equations:
$$
\nabla \cdot \sigma = 0
$$
where $\sigma$ is stress (contains both the mechanical and thermal stresses).

I currently am using the finite volume method to solve the above equation, and the integral form is
$$
\int \int (\sigma \cdot \hat{n}) dA = 0 \ \ \ \ \text{   applied divergence theorem}\\
\int \int (\sigma \cdot \hat{n}) dA \approx \sum_{f=1}^N \sigma_fA_f\hat{n}_f
$$
where the integral is approximated using the midpoint rule that reads as the sum of the product of the face stress, face area, and face normal over all faces of a cell. The face stress is computed by taking an average of the neighboring cell's stress.

The stress uses the duhamel-neumann relationship for an anisotropic material (note  that I am currently testing the code on a simple isotropic material)
$$
\sigma = C:(\epsilon_m - \epsilon_t)
$$
Where $\epsilon_m$ is the mechanical strain that is computed using the linear strain displacement relationship and $\epsilon_t$ is the thermal strain, computed using
$$
\epsilon_t = \alpha (T-T_0)
$$
where $\alpha$ is the CTE.

I can formulate the above system of equations as a linear system $Au=b$, where $u$ are the displacements to be solved for. $b$ is the forcing and essentially all the thermal stress term is contained inside $b$ (thermal stress has no impact on $A$).

The good news is that if I assign the same thermal stress to the entire domain (i.e., same temperature at all the cells), then I get a solution that is "correct." By correct, I mean it matches commercial software (ANSYS). So something is correct.

However, if my temperature distribution across the domain is not the same, then I get solutions that not only do not come close to matching ANSYS, but also do not match my intuition (certain cells with larger temperature changes are expanding less or even contracting(!), as indicated by the computed displacement field, than parts with smaller or no temperature changes.

I have checked my code and it seems I am implementing everything correctly, so I am wondering if my numerics are incorrect. But I am getting the correct solution if the whole field is initialized to the same temperature, so something must be correct. e.g., does adding the thermal stress term somehow affect how I should approach the discretization? I would think not. This term is easy and it's really just a forcing term..

As previously mentioned, the thermal stress is manifested mathematically in the forcing vector $b$ and has no impact on $A$. So changing the temperature just affects the magnitude of my forcing vector, and my $A$ remains the same. So this tells me that my $A$ matrix is formulated correctly, or that it could be formulated incorrectly, but in a way that has not affected my solutions for other problems. I think the latter is very unlikely because I have also tested my code on cantilever and various other bending/compression problems (with NO thermal stress), and I get solutions that match well with analytical solutions or with commercial software. I have also done some manufactured solutions studies to look at convergence and the tests succeeded.

It's just this thermal stress term that's causing me issues, and only when the temperature is not uniform across the discretized domain.

Any suggestions? I'm happy to provide more info if needed.

Edit #2:
Here is the particular geometry I am working with right now:

The length is in the x direction. I apply a uniform temperature of $1900K$ to all the cells. The following picture is the displacement contours for the y-direction, which is the more refined direction in the above picture.

The max y-displacement from the contours for the +y and -y surfaces are about 0.0006m and -0.0006m. This material is steel, using CTE =  $1.2\times 10^{-5} K^{-1}$. The dimension of the block is about 0.6098m by 0.0047625m by 0.00508m (sorry multiple decimals. The original grid was created in inches) in the x,y,z directions. The mesh step size in the 3 directions is constant, at 0.03048m (1.2inches) by 0.004318m (0.17 inches) by 0.0254m (1 inch).

So here I get a solution that matches what ANSYS gives me.

Next I proceeded to run a different simulation where I only specify a temperature of 1900K in the domain in the top two layers (see the region in the white boxes below) and the rest of the domain is at 298K.

In this simulation where only part of the domain experiences a temperature change, I get a solution that is way off from expected or from what ANSYS gives me.

Edit #3: Adding displacement contours for the aforementioned test case.

Here are the 3 displacement contours for the x-,y-, and z-displacements for the second case. My intuition tells me that the x- and z- displacement contours show trends that are expected (not sure if the magnitudes are correct, but at least the trends are reasonable). But I would not expect the y-displacement to behave the way it is shown in the y-displacement contours.


问:
Element Preconditioner
Im just working on a preconditioner for the linear equation system $Ax = b$ arising in FEM for elliptic PDE. $A$ is a s.p.d Matrix with real valued entries. I read something about the element by element preconditioner and now I have a question about the following:

Defining a preconditioner $P$ such that $P \approx A$. One can solve the equation system (left preconditioning) instead:

\begin{equation}
P^{-1}Ax = P^{-1}b
\end{equation}

I read that some people choosing the preconditioner as the Assembly of Element Preconditioners.

\begin{equation}
P =\prod\limits_{e =1}^n P_e
\end{equation}

where $P_e$ is defined for every element $e$ as
\begin{equation}
P_e = I + diag(A_e)^{-1/2}\cdot(A_e - diag(A_e))\cdot diag(A_e)^{-1/2}
\end{equation}

Now to calculate $P^{-1}$ you just need to assembly all $P_e^{-1}$.

My question is now: If I use a Cholesky decomposition would for $P$ such that $P = C\cdot C^T$ would it be the same as if I use the decomposition for every element and assemble afterwards:??
\begin{equation}
P_e = C_e C_e^T \rightarrow C\cdot C^T = \prod\limits_{e = 1}^n C_e C_e^T
\end{equation}

I dont have the feeling that this holds because of overlapping Elements. But Im not sure.

By the way my goal is to use a decomposition to solve

\begin{align}
C^{-1}AC^{-T}y &= C^{-1}b \\
x &= C^{-T}y
\end{align}

Title: Element Preconditioner


问:
Why is the method of im2col with GEMM is more efficient than the method of direction implementation with SIMD in CNN
The convolutional layers are most computationally intense parts of Convolutional neural networks (CNNs).Currently the common approach to impement convolutional layers is to expand the image into a column matrix(im2col) and perform and perform Multiple Channel Multiple Kernel (MCMK) convolution using an existing parallel General Matrix Multiplication (GEMM) library. However im2col operation need load and store the image data, and also need another memory block to hold the intermediate data.

If I need to optimize the convolutional implementation, I may choose to direct implementation with SIMD instructions. Such method will not incur any memory operation overhead.

the benefits from the very regular patterns of memory access outweigh the wasteful storage costs.

From the following link, at the end of the link

https://petewarden.com/2015/04/20/why-gemm-is-at-the-heart-of-deep-learning/

So I hope to know the reason. May floating point operations require more instruction cycle? or the input image is not much large, so it may residue in the cache and the memory operations don't need access DDR and consume less cycles.


问:
Factorize laplacian in terms of first derivative matrix
I am trying to factorize the following Laplacian matrix in terms of $ D^TD$, D is the first derivative matrix.
The tridiagonal form of the secon derivative matrix using Neumann boundary condition is given by,

$\begin{bmatrix} -1 & 1 & 0  & 0 & 0 \\ 1 & -2 & 1 & 0 &0\\ 0 & 1 & -2 &1&0\\0&0&1&-2&1\\0&0&0&1&-1\end{bmatrix}$.

In order to write the above in terms of the gradient( first derivative) operator,  I tried factorizing the Laplacian using Cholesky factorization. But, I couldn't succeed in factorizing.

Any suggestions on how to find a decomposition of the form $D^TD$?


问:
Numerically solving a system of parabolic PDEs and 1st order ODEs
I'm trying to solve the following system of differential equations numerically. What are the available finite difference approaches and matlab solvers to solve such a system? Other approaches to solve the problem are welcome as well. The first two equations are quasi-steady and the rest of the equations have time derivatives in them. I'm trying using matlab's pdepe solver, but since the boundary conditions have to be specified even for the ODEs which have time-derivatives, I didn't think it as a viable option.

$$0= \rho \frac{\partial}{\partial z}\left(D \frac{\partial c_1}{\partial z} \right) - \it{source_1}\, $$

$$ 0= \rho \frac{\partial}{\partial z}\left(D \frac{\partial c_2}{\partial z} \right) - \it{source_2}\,$$

$$\frac{\partial c_3}{\partial t} = source_3$$

$$\frac{\partial c_4}{\partial t} = source_4$$

$$D = A_1 c_1  + A_2 c_2 + \dots + A_n c_n \,,
\rho = 1 (say)\,.$$

The source terms for any particular equation (or concentration) contain dependence on the other concentrations.

OwnerDisplayName: mnatch


问:
Fast calculation of $A^T B$
I need to compute a matrix-matrix product, $A^T B$, where $A$ is $n \times r$ sparse, and $B$ is $n \times q$ dense. The number of rows $n$ is far larger than both $r$ and $q$. In fact $n$ is so large I cannot store the entire matrix $B$ in memory, so I build it one row at a time and update the matrix product $A^T B$ accordingly. The way I currently do this is to store $A$ in CSR format and then apply the following algorithm:

```<code>loop i = 1, ..., n     // loop over rows of B
  bi = build_row(B, i) // construct row i of B
  loop p = A->row_ptr[i], ..., A->row_ptr[i+1] - 1 // loop over non-zero entries in row i of A
    DAXPY(A->data[p], bi, ATB(A->col_ind[p], :))
</code>```

This algorithm seems to have complexity $O(nkq)$, where $k$ is the number of non-zero entries in a row of $A$ (this is the same for all rows of $A$). This algorithm works fine, but is quite slow so I am looking for ideas on how to speed it up.

Dense BLAS algorithms often take advantage of blocks of data instead of operating on one row at a time, so my thinking was to do the same (instead of building B one row at a time, build a block of rows, where the block size is still small enough to fit in memory, and then apply a matrix multiplication algorithm). For this type of method, I don't think CSR format would work, since we would need the columns of $A$ in sequential memory locations, so CSC format might work better. A possible algorithm (with A in CSC format) is:

```<code>loop i = 1, ..., num_blocks
  row0 = (i-1) * block_size           // row of start of block in B
  row1 = row0 + block_size            // row of end of block in B
  blockB = build_block(B, row0, row1) // build B(row0:row1,:)
  loop j = 1, ..., r // loop over rows of A^T
    loop p = A->col_ptr[j], ..., A->col_ptr[j+1]-1 // loop over non-zero entries in row j of A^T
      if (row0 <= A->row_ind[p]) AND (A->row_ind[p] <= row1) // if this row of A is in the right range, perform DAXPY operation
        DAXPY(A->data[p], blockB(A->row_ind[p]-row0, :), ATB(j, :))
</code>```

I think the above algorithm has complexity $O(k r q \times NB)$, where $k$ is the number of non-zero entries in a column of $A$, and $NB$ is the number of blocks (num_blocks).

I haven't yet coded up this method, since I'm not sure it will actually be faster. Does anyone have any suggestions on how to make a fast algorithm to do this? Perhaps there is another sparse matrix storage format which would be better to use?

Title: Fast calculation of $A^T B$


问:
Approximate the largest simplex in N-dimensional Delaunay triangulation
I am working on determining the spatial information of a set of $M$ points in $N$-dimensional space. It is well-known that the construction of Delaunay triangulation is expensive in high dimensional space. Thus, I would like to only find the largest simplex of $N$-dimensional Delaunay triangulation (DT). I am going to sample a new point inside that simplex by each iteration, which means that the largest simplex would become smaller and smaller by iteration.

Is there any idea or references to computing the largest simplex in DT? It may not be exactly the largest simplex after computing DT. But it should at least (1) approximate the largest simplex within a tolerance; (2) the largest simplex that this method computed should become smaller and smaller by each iteration.


问:
Why the MIRACLE of Lanczos/CG-like?
Lanczos/Arnoldi/Rietz/CG-like algorithm share the same core strategy...
In each, a little miracle appears, most of the Gram-Schmidt inner products are zeroes! In others words, new direction need only get orthonormalized with the one previous updated search direction (or two in Lanczos case), no need to check other previous directions (because of zeros). Rather a big miracle when $A$ is sparse and $n = 10^6$.

In case of non-linear inverse/eigenproblem, this property is not true, but depending of the "nonlinearity rate", the values could be small and easily neglected or threshold since the iterative nature of the algorithm assumed that the residual must decrease in each step. So this property remains useful.

Any suggestions are welcome for understanding why the new direction is "already" orthogonal to the set of all (previous+1) directions.
Math is ok, but there is a deep principle behind this thing.
It looks like a hidden tridiagonalization process that happens at the same time. Is there an easy way to handle the trick?

Title: Why the MIRACLE of Lanczos/CG-like?


问:
Minimum axis aliged bounding box of convex polytope
I need to compute a $n-$dimensional integral with $n<10$ on a convex polytope. Since most numerical integration libraries (e.g. Cuba) expect the function to be integrated defined inside an axis aligned bounding box (AABB), I need to compute a (minimum volume) AABB enclosing the polytope. Since the polytope is defined by a set of linear inequalities (H-polytope), the easiest way to compute the AABB would be converting the polytope to the V-form using the Fourier–Motzkin elimination (FME) algorithm. The problem is that the number of vertices generated by FME grows so quickly to make the method not practical. So my question is: is there a way to compute the AABB directly from the H-polytope while avoiding the conversion to V-polytope?


问:
Why do stabilized formulations for the Navier-Stokes equation maintain the convergence rate for high order polynomial interpolation?
I have a quick questions which has been troubling me lately.
When reading the FENICS Finite Element Book they assess various approaches to solver the Stokes equation. Obviously, they discuss the notion of elements that satisfy the LBB conditions, but they also discuss stabilized method. They briefly explain that the core of stabilization method is to add a term to the pressure equation which alters the saddle point problem at the consequence of making the velocity "less divergence free" in layman's term. This is discussed at page 397 of the pdf.

The resulting weak form reads :

They noticeably insist on the fact that consistency is lost for orders q>1 which makes sense because $\delta$ is of order $h^2$, consequently it would appear to me that the convergence rate for pressure becomes bounded by this term if the order of the pressure interpolation is higher than 2.

Later on in the book they look at the order of convergence of the various scheme as well as the mass loss on a manufactured solution. They look at various polynomial order for the various elements and show the convergence order. For the stablized method, they observed that consistency is maintained for the L2 error on the velocity, for instance using 4th order Lagrange polynomial interpolation, they obtain a order of convergence of 4+1.
This is shown in the following table (see the STAB column):

I do not understand how this occurs, considering that the term that is added to the continuity equation is by itself of order 2 in L2 space (due to it being $h^2$). What is wrong in my reasoning? I am obviously missing some details and something is wrong in my understanding of stabilization + high order interpolation. Is consistency maintained for velocity but lost for pressure (that is the order of convergence for pressure is reduce to second order, but not for velocity)?
Thank you very much!


问:
Solving multiple linear programs with same constraints but different objective
I have ~30 non-negative variables and 24 equations and I want to find out the upper and lower bound for each variable. Feasible solutions are guaranteed.

So for each variable, I solve two LP problem, one for min and one for max

$$
\begin{array}{rl}
\min & c^T x\\
\textrm{s.t.} & Ax = b\\
&x\geq0
\end{array}
$$

where $c=(0, 0, \dots, 0, 1, 0, \dots, 0)^T$ and $(0, 0, \dots, 0, -1, 0, \dots, 0)^T$.

So in total I have 60 LPs with exactly the same constraint, but slightly different objective. I have to do the same thing for each day's data for a few thousand days.

Is there a way to somehow reuse the intermediate results to speed up the calculation? Is it possible to implement using existing Python or R package?

Edit: Just to add that for each day, only the rhs $b$ is different, everything else stays the same.


问:
Guide for finite-difference schemes for Hamilton-Jacobi-Bellman Equations
I need to solve a simple, low-dimensional Hamilton-Jacobi-Bellman equation.

Is there a simple guide for doing this numerically using finite-difference schemes? I found a few research articles discussing higher order methods or discussing theoretical properties of numerical schemes but nothing that really helped me, as a newcomer without academic aspiration, to implement an actual algorithm. I really don't care about runtime and am perfectly fine with something that takes a whole day for 5% relative accuracy.

When I tried a simple forward Euler method, I got bizzare oscillations even though I tried to evaluate derivatives upstream. What upstream is depends on the solution of the Hamiltonian, of course, which itself requires evaluation of the derivatives; hence, there is some circularity here and I guess I got it wrong. Alternatively, I might have messed up the boundary conditions.

As an example of what I need to solve:

$$
u\colon [0,T]\times \mathbb{R}^d\to \mathbb{R}
$$
$$
u_t-H(u,\nabla u)=0
$$

$$
H(u,\nabla u)=\min_{\alpha_i\geq 0\\ a\cdot \alpha=1} b\cdot \alpha+(c-\alpha)\cdot \nabla u
$$
$$
u(0,\cdot)=0
$$

I tried the following iteration over finite two-dimensional arrays $(u_n, \nabla u_n)$, $0\leq n\leq N$ (with $\Delta t:= T/N$):

(For each array index) solve the Hamiltonian minimization problem to get $H(u_{n},\nabla u_{n})$ and $\alpha$.

(For each array index) define $u_{n+1}:=u_n+\Delta t H(u_{n},\nabla u_n)$

(For each array index) define $\nabla u_{n+1}$ by differencing $u_{n+1}$ in upstream direction, as defined by $\alpha$

At the boundary, I used $u_{x_ix_i}:=0$, i.e., I whenever the upstream method required me to access a non-existent entry of $u$, I evaluated downstream instead.


问:
Is it possible to obtain a 'relaxing lengths' for 8 node hexaedral element
I am using 3D FEM with 8-node brick elements to model a certain type of growth/expansion, which is not a plastic deformation. First, I apply a pressure/load and I get displacement and the new position of nodes and calculate stress/strain. Next, I would like to change the 'edges' relaxing length (different edges by different factors) of an 8-node hexahedral element. Such change in relaxing length would be due to addition of a new material, resulting in expansion of a mesh.

My question is, whether in such a model I can calculate some sort of relaxing length for a brick element edges after application of pressure and consequently change it to get new nodal positions. I can imagine doing so for a 1D string, but I am not sure if it makes any sense for a 3D hexahedral element.


问:
Improving convergence of Jacobi iteration to Schur form
I'm using SIMD processor arrays to compute the eigen-decomposition for large numbers of small (up to $32\times 32$) matrices. For assorted technical reasons, Jacobi iteration maps well to the SIMD architecture and gives excellent results for diagonalisation of symmetric, Hermitian and normal matrices using unitary similarity transformations.

I've extended the technique to compute the Schur triangular form of general complex matrices using the same unitary similarity transformations in order to obtain the eigenvalues from the diagonal. Again, convergence to the Schur form is typically excellent for most random complex matrices up to $32\times 32$ size but I observe very poor convergence for a small number of random matrices. As I increase the size of the number of random matrices tested to thousands, tens of thousands, millions, the number of Jacobi iterations required increases indefinitely.

Is there some form of pre-conditioning that can improve convergence? I've tried an initial reduction to upper Hessenberg (i.e. almost triangular) with no improvement.


问:
Are there known accuracy issues between 2D axisymmetric and 3D solutions?
In my full 3D solutions I am solving for the potential throughout a $100\times 200\times 200$ grid. Inside is a ring electrode set to -5V via a Dirichlet boundary condition, and surrounded on all sides by Dirichlet boundary conditions at 0V.

When I solved for this setup in full 3D vs solving for the same setup in 2D axisymmetric form, I get different results. In the images below, one can denote the darker blue color, in the very center of each of the two solutions.

For the 2D axisymmetric setup, the left, top, and right all have the same 0V Dirichlet boundaries as the 3D setup except for the Neumann boundary on the bottom where the axis of symmetry for the problem resides.

Below you can see two images. In both images, the object on the left side is a 'slice' of the 3D solution, and the object(s) on the right is the solution for the 2D axisymmetric solution. The first image shows the 2D axisymmetric solution, whereas the second image has a 'mirror' image of that solution placed just below its counterpart, to make comparing the two solutions (3D vs 2D axisymmetric) that much easier to see.

So, is it in any way common for such a marked difference in the results of 3D vs 2D axisymmetric results for elliptic PDEs?

UPDATE: I would like to thank Bill and Wolfgang for their constructive questions.

Knowing that the 2D Axisymmetric and 3D solutions should be the same, I modified my 3D setup and used a $100\times199\times199$ grid instead, with $(i, 99, 99)$ as the axis for the ring electrode, and setting Dirichlet boundaries at a radius of $100$ from said axis. The solutions in 3D now match my 2D axisymmetric solution.

Below are the images showing the now matching solutions.


问:
2d wave equation with finite differences blowing up
I am (naively) trying to solve the 2d wave equation with finite differences. But the system blows up instantly.

For simplicity I set the constant $c=1$, then I am left with
$$\Delta u =u_{tt}.$$

I set up the FDM with

```<code>N = 9
u = np.zeros((N**2,1))
u.reshape(N,N)[N//2,N//2] = 1
</code>```

Now, at the time t=0 the function looks like this

```<code>u_0 =
[[0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 1 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]
 [0 0 0 0 0 0 0 0 0]]
</code>```

To which I apply the five point diffence stencil Lh

```<code>import scipy.sparse as sp
import numpy as np
h = 1/(N+1)
T = np.diag(N*[-2]) + np.diag((N-1)*[1],k=-1) +np.diag((N-1)*[1],k=1)
spt = sp.coo_matrix(T)
spi = sp.eye(N)
Lh = (sp.kron(spi,spt) + sp.kron(spt,spi) ) / h**2
</code>```

Since Lh is my discretized Laplacian, for every timestep I add the acceleration that it gives me to some book-keeping vector v to change the velocity accordingly:

```<code>v = np.zeros_like(u)
print(u.reshape(N,N).astype(int))
for i in range(3):
    v += Lh.dot(u)
    u += v
    print(u.reshape(N,N).astype(int))
</code>```

Here is what happens:

```<code>[[   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0   99    0    0    0    0]
 [   0    0    0   99 -398   99    0    0    0]
 [   0    0    0    0   99    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]
 [   0    0    0    0    0    0    0    0    0]]
</code>```

In the next timestep:

```<code>[[     0      0      0      0      0      0      0      0      0]
 [     0      0      0      0      0      0      0      0      0]
 [     0      0      0      0   9999      0      0      0      0]
 [     0      0      0  19999 -79699  19999      0      0      0]
 [     0      0   9999 -79699 198800 -79699   9999      0      0]
 [     0      0      0  19999 -79699  19999      0      0      0]
 [     0      0      0      0   9999      0      0      0      0]
 [     0      0      0      0      0      0      0      0      0]
 [     0      0      0      0      0      0      0      0      0]]
</code>```

What is the problem here? I have not set up a boundary condition to my solution. Maybe I should?

greetings


问:
Where does the seemingly official number of certain algorithms come from?
There are a lot of algorithms which seem to have been supplied an official number, such as Algorithm 76, Hierarchical clustering using the minimum spanning tree. Another example is Algorithm 123, Single Transferable Vote by Meek's Method. These numbers seem to refer to some official registry, but I haven't been able to find out what this registry is or how numbers are assigned. It does not seem to be related to the journal publishing the algorithm. Does anyone know how this system works and where one can read more about it?</p>


问:
Gmsh for 3D volume with inclusions
In an attempt to create three-dimensional volumes with inclusions in Gmsh I stumble upon a problem which was non-existent in the two-dimensional case.

I'm using the OpenCASCADE geometry kernel because this facilitates convenient boolean operations to remove parts of the inclusions outside my considered domain box.

In a nutshell the problem boils down to:

whenever I incorporate inclusions I perform a unite step, such that the different geometrical features are recognized but the interfaces
are joined in the mesh. This works like a charm in 2D, however, when
I do the equivalent in 3D, my volumetric geometrical information is
lost and volumes are actually united completely removing the interface.

I can graphically depict this by performing the 2D mesh step, in both 2D and 3D:

What the above showcases is the intended behaviour in 2D where (continuous) interfaces are retained, but although in 3D the sphere locations are still visible from the outside, the internal surface is no longer existent, i.e. the sphere outside is not meshed, only the cut section.
I should note that I believe I need the unite step to tie the different geometrical features together, that I've tried many different minor variations to produce 3D results like the one depicted above but none resulted in the intended mesh and that my web search was unhelpful.

I hope someone has a concrete idea or even a hunch on where to find the solution for the issue.

For completeness I'm adding the simplified .geo examples used to obtain the above graphs.

For 2D:

```<code>SetFactory("OpenCASCADE");

Point(1) = {0, 0, 0};
Point(2) = {2, 0, 0};
Point(3) = {2, 2, 0};
Point(4) = {0, 2, 0};

Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};

Line Loop(1) = {1, 2, 3, 4};

Plane Surface(1) = {1};

Periodic Line{1} = {-3};
Periodic Line{2} = {-4};

Disk(2) = {0, 0, 0, 0.5};
Disk(3) = {2, 0, 0, 0.5};
Disk(4) = {0, 2, 0, 0.5};
Disk(5) = {2, 2, 0, 0.5};
Disk(6) = {1, 1, 0, 0.5};

f1() = BooleanDifference{ Surface{1}; }{ Surface{2:6}; };
f2() = BooleanIntersection{ Surface{1}; Delete; }{ Surface{2:6}; Delete; };
BooleanUnion{ Surface{ f1() }; Delete; }{ Surface{ f2() }; Delete; }
</code>```

For 3D:

```<code>SetFactory("OpenCASCADE");

incRadi = 0.5;
boxSize = 1.5;

Box(1) = {0,0,0, boxSize,boxSize,boxSize};

Periodic Surface{2} = {1} Translate{boxSize,0,0};
Periodic Surface{4} = {3} Translate{0,boxSize,0};
Periodic Surface{6} = {5} Translate{0,0,boxSize};

Sphere(2) = {0, 0.5*boxSize, 0.5*boxSize, incRadi, -Pi/2, Pi/2, 2*Pi};
Sphere(3) = {boxSize, 0.5*boxSize, 0.5*boxSize, incRadi, -Pi/2, Pi/2, 2*Pi};

f1() = BooleanDifference{ Volume{1}; }{ Volume{2:3}; };
f2() = BooleanIntersection{ Volume{1}; Delete; }{ Volume{2:3}; Delete; };
BooleanUnion{ Volume{ f1() }; Delete; }{ Volume{ f2() }; Delete; }
</code>```

Title: Gmsh for 3D volume with inclusions


问:
Identify the components of the (weak form) PDE in structural mechanics
I am trying to identify the weak form of PDE in structural mechanics. I read a lot of papers where they are using the elliptic boundary value problem

\begin{equation}
 \int\limits_{\Omega} \delta \epsilon : \sigma \,d\Omega - \Big[\int\limits_{\Omega} \delta u \cdot b \, d\Omega + \int\limits_{\Gamma_N} \delta u \cdot t \, d\Gamma\Big] = 0
\end{equation}
where $b$ represent the body forces, $t$ is the traction vector acting on the Neumann boundary, stress $\sigma$ and $u$ the displacement. With $\delta \epsilon$ and $\delta u$ they define virtual strain/displacement.

To be honest, I am not used to mechanics. For example I am using the general elliptic PDE:

\begin{align}\left\{\begin{aligned}
div(A\nabla u) = & \,f && \text{ in } \Omega   \\
u  = & \,g_D && \text{ on }\Gamma_D \\
A \nabla_n u = &\, g_N && \text{ on }\Gamma_N \\
  \end{aligned}\right.
\end{align}
with $A$ a bounded self-adjoint and elliptic linear operator. The weak form this PDE reads: find $u\in H^1_{\Gamma_D}(\Omega)$ for given functions $f\in L^2(\Omega)$ and $g_N \in L^2(\Gamma_N)$ such that

\begin{equation}
\int\limits_{\Omega} \nabla u A \nabla v  \, d\Omega = \int\limits_{\Omega} v^T f \, d\Omega + \int\limits_{\Gamma_N} v^T g_N \,d\Gamma \qquad \forall v \in H^1_{\Gamma_D}(\Omega).
\end{equation}

My goal now is to identify the weak form in structural mechanics such that my general weak form will be the same. As I can see, first I need to define $v^T = \delta u$, $t = g_N$ and $b = f$. But I am not sure about the other terms. Is it correct that $\delta \epsilon = \nabla v = \nabla \delta u$ and $\nabla u = \sigma$? I cannot see the connection between these definitions. And how do I use the elliptic operator $A$ in this context? Can someone bring light into the darkness?

My new Idea is: Use that strain is equal to the derivatives of the displacement such that $\epsilon = \nabla u$. Next is to define the virtual strain as the derivative of the virtual displacement such that $\nabla v = \nabla \delta u = \delta \epsilon$. By Hookes law it is now:

\begin{equation}
\nabla v : A \nabla u = \nabla \delta u : A \nabla u = \delta \epsilon : A \epsilon = \delta \epsilon : \sigma
\end{equation}


问:
Comparing sum of floating points
I am currently working on a numerical algorithm involving a lot of floating point arithmetic, involving some badly conditioned problem sets.

I am using the relation $|x - y| / (\max(|x|, |y|, 1)) \leq \epsilon$. Unfortunately, the corresponding relation $=\sim$ is not invariant with respect to additions, i.e. $x =\sim y$ does not imply $x - y =\sim 0$.

As a counter example I have $x = 1e+20$, $y = 1.0000000000000082e+20$, and therefore $x - y = -819200$ (wrt. floating point arithmetic). I think that this is a pretty typical case of cancellation...

My questions are:

Is there a suitable relation?

How should I compare $x + y =\sim z$? I am using the statement
in an assertion which consequently fails on certain instances.
Is it advisable to use ($x + y =\sim z$ or $x =\sim z - y$)
instead?

Edit:

What am I trying to achieve: I am implementing a solver for nonlinear programming problems, and running the implementation on the cutest suite of test problems. In the specific step, I am asserting that

$$ c(x) + J(x)d = \sim 0, $$

where $c(x)$ is a vector of constraint values, and $J(x)$ is the constraint Jacobian. The problems in the test set vary widely in terms of size and type. Therefore, I'm afraid that there is not much specific knowledge.

What is more, the test suite is known to be badly scaled, i.e. containing constraints of the type $1000*x^2 + 0.001 * y^2$. Therefore, the Jacobian is likely ill-conditioned.

Title: Comparing sum of floating points


问:
Algebraic recursion of Hermite polynomials in SymPy
I want to obtain the algebraic values of, for example, Hermite polynomials using SciPy but in a recursive manner.

Using Maple, for example, these can be defined as

```<code>H(0,x):=1;
H(1,x):=2*x;
for n from 1 to 6 do
    H(n_1,x):= 2*x*H(n,x) - 2*n*H(n-1,x);
    print( simplify(%) );
end do:
</code>```

I can find no way in SciyPy to do this algebraically. Please can you help?


问:
Solve system of polynomial equations with Python
I have 5 at most 4th order polynomials in 5 variables,

$$p_i(x_1,x_2,x_3,x_4,x_5) \qquad i = 1, \ldots, 5$$

where all coefficients are either rational or floating point. I'd would like to get the roots of the simultaneous equations $p_1 = 0, p_2 = 0, p_3 = 0, p_4 = 0, p_5 = 0$ and I'm partially able to do this quite easily with non-linear solvers using python scipy. Of course these non-linear solvers don't give me all roots just a couple of them if I choose the starting points carefully. In any case there is no guarantee that I'm getting all of them. But since we are talking about polynomials there are known algorithms for reducing them to polynomials of a single variable and the full set of roots of those can be found easily (see Wiki).

What I'm wondering is if there was any of this implemented in Python/SciPy. What I've found so far is finding the roots of a polynomial (single variable) and general non-linear solvers.


问:
Nonlinear least squares and regularization
Consider the nonlinear least-squares minimization of a vector of $n$ residuals $\mathbf{f}$ in $p$ parameters $\mathbf{x}$:
$$
\min_{\mathbf{x}} || \mathbf{f}(\mathbf{x}) ||^2
$$
This can be done with a standard Gauss-Newton or Levenberg-Marquardt approach, using the $n \times p$ Jacobian matrix
$$
J_{ij} = \frac{\partial f_i}{\partial x_j}
$$
Now suppose I want to regularize the solution, so I want to solve
$$
\min_{\mathbf{x}} || \mathbf{f}(\mathbf{x}) ||^2 + \lambda^2 x^T C x
$$
for some symmetric $p \times p$ matrix $C$. It seems to me there are two ways to proceed.

Option 1

If $C$ happens to be SPD, then we can do a Cholesky decomposition $C = U^T U$ where $U$ is upper triangular. Then we get
$$
\min_{\mathbf{x}} || \mathbf{f}(\mathbf{x}) ||^2 + \lambda^2 || U \mathbf{x}||^2
$$
Here we can form the augmented residual vector and Jacobian:
\begin{align}
\tilde{\mathbf{f}} &= \pmatrix{ \mathbf{f} \\ \lambda U \mathbf{x}} \\
\tilde{J} &= \pmatrix{ J \\ \lambda U}
\end{align}
and the problem becomes
$$
\min_{\mathbf{x}} || \tilde{\mathbf{f}}(\mathbf{x}) ||^2
$$
for which we can apply a standard Gauss-Newton or LM code to solve.
The augmented system will have $n + p$ residuals and a $(n+p) \times p$ Jacobian matrix. For a large number of parameters $p$ (say several thousand), this could add significant cost to the nonlinear iteration.

Option 2

Alternatively, suppose $C$ is symmetric but not positive definite, then we can't perform a Cholesky factorization. In this case we can still proceed by forming the augmented system:
\begin{align}
\tilde{\mathbf{f}} &= \pmatrix{ \mathbf{f} \\ \lambda \sqrt{x^T C x}} \\
\tilde{J} &= \pmatrix{ J \\ \lambda \frac{x^T C}{x^T C x}}
\end{align}
In this case, we now have $n + 1$ residuals and the Jacobian is $(n+1) \times p$.

My question is whether Option 2 is superior to Option 1? Option 2 requires factoring a smaller Jacobian matrix which could result in significant speed improvements when $p$ is very large. However I wonder if Option 1 might converge faster to the solution, since more "information" is stored in the residual vector and Jacobian. Does anyone have any insights into this?

For a long time I used Option 1, but recently I encountered a problem where $C$ is symmetric but not positive definite, so I am forced to consider Option 2. Just wondering if Option 2 might have worse convergence properties.

EDIT: it seems there is a third option which may or may not be useful:

Option 3

The minimization problem can be put in the form:

$$
\min_{\mathbf{x}} \tilde{\mathbf{f}}^T W \tilde{\mathbf{f}}
$$
where
$$
W = \pmatrix{ I & 0 \\ 0 & C}
$$
and
$$
\tilde{\mathbf{f}} = \pmatrix{ \mathbf{f} \\ \lambda x}
$$
In this form, it may not be possible to use an off-the-shelf LM software due to the $W$ matrix, but certainly a specialized algorithm could be developed to solve this problem. Again I am not clear if this offers any advantage over Option 1 or 2.


问:
Any way to avoid catastrophic cancellation when computing the discriminant of a quadratic function?
Homework disclaimer...

The task:

We are using the following algorithm to solve the quadratic equation $x^2+2px+q=0$:

$x_1=|p|+\sqrt{p^2-q}\mathtt{;}$

$\mathtt{if}\,p>0\,\mathtt{then}\,x_1=-x_1\mathtt{;}$

$x_2=q/x_1\mathtt{;}$

Point out two places where running this algorithm can meet significant problems. How to mitigate these problems?

That's my guess: The first place is this: $p^2-q$, we're risking catastrophic cancellation here; and the second place is this: $q/x_1$, we're getting a NaN if $x_1$ is $0$.

The second one seems kind of easy: the only way for the absolutey larger root to be $0$ is when both roots are $0$, which entails so $p=q=0$; we can test for this case.

I don't know how to mitigate the first problem, though. I can't think of a way to solve a quadratic equation without computing its discriminant, and computing its discriminant does require subtraction. Which, of course, risks catastrophic cancellation.

Is there any way to avoid catastrophic cancellation when computing the discriminant of a quadratic function?

EDIT: As per request from comments, moving additional info from answer to question:

Hah. Regarding $x_1=|p|+\sqrt{p^2-q}$. Another student was arguing today that this cannot be an issue: According to him, catastrophic cancellation can only occur if $p^2\approx q$; but then $\sqrt{p^2-q}\approx 0$ and moreover $|p|\gg\sqrt{p^2-q}$; so the discriminant itself and any errors introduced by computing it loose significance when computing $x_1$. The Professor said he would have to think about it; I'll post an update here if/when he's done thinking.

However, there is one more problem I failed to notice. Namely, if $p$ is large, then $p^2$ may overflow; but then $\sqrt{p^2-q}$ will also overflow even if its accurate value would fit in the range of the number format; so $x_1$ will be spuriously set to $\infty$. To avoid this, the Professor said, $p$ needs to be rescaled if its large enough.


问:
Is there a numerically stable way to take $\epsilon \rightarrow 0$ in integrals of the form $\int \frac{f(x)dx}{x+i\epsilon}$?
The Sokhotski-Plemelj theorem states,
$$\lim_{\epsilon\rightarrow 0^+}\int_a^b\frac{f(x)dx}{x+i\epsilon} = \mathcal P \int_a^b \frac{f(x)dx}{x} - i\pi f(0). $$

Is there a numerically stable way to take this limit, without explicitly using the above theorem? I find that if I use $f(x)=1$, $a=-b$, I almost converge to $-i\pi$, but not before things become unstable.


问:
Solution of constrained system of ODEs
Can someone point me in a direction to solve this kind of integral constrained system of ODEs.

\begin{align}
&\int_0^{1/2}\dot{y}^2(t)=p\\
&2\lambda_1\ddot{y}(t)+\pi cos(\pi y(t))=0\\
&y(0)=0,y(1/2)=1/2
\end{align}

I have reduced it to 1st order:
\begin{align}
&\int_0^{1/2}x^2(t)=p\\
&\dot{y}=x \\
&2\lambda_1\dot{x}(t)+\pi cos(\pi y(t))=0\\
&y(0)=0,y(1/2)=1/2
\end{align}

but its still not suitable for an ODE solver. Any help will be appreciated.

Edit:$p$ is a known constant and $\lambda_1$ is an unknown constant.

Disclaimer: This is a cross-post and has some good answers here on MathSE.

Title: Solution of constrained system of ODEs


问:
Nonlinear system with diagonal nonlinearity
Consider a nonlinear system of the form $\boldsymbol{f}(\boldsymbol{x}) = \boldsymbol{0}_{\mathbb{R}^n}$ for $\boldsymbol{x} \in \mathbb{R}^n$, where the function $\boldsymbol{f}$ is given by
\begin{equation}
\boldsymbol{f}(\boldsymbol{x}) := \boldsymbol{g}(\boldsymbol{x}) + \boldsymbol{\underline{A}} \boldsymbol{x} + \boldsymbol{b} \in \mathbb{R}^n, \quad \boldsymbol{x} \in \mathbb{R}^n.
\end{equation}
Here, the function $\boldsymbol{g}$ is differentiable. Moreover, the functions $g_i$ depend only on $x_i$, i. e. $\frac{\partial g_i}{\partial x_j} \equiv 0$, $i \neq j$. The matrix $\boldsymbol{\underline{A}} \in \mathbb{R}^{n \times n}$ is large, sparse, non-symmetric, and invertible.

Can we take advantage of the purely diagonal nonlinearity in $\boldsymbol{f}$ somehow when solving the nonlinear system?

The assembly of the Jacobian matrix $\boldsymbol{\underline{J}_f}(\boldsymbol{x}) = \boldsymbol{\underline{J}_g}(\boldsymbol{x}) + \boldsymbol{\underline{A}}$ is cheap since $\boldsymbol{\underline{J}_g}(\boldsymbol{x})$ is a diagonal matrix. Therefore, Newton's method appears to be a natural choice, but other methods are also acceptable.

If we use Newton's method, then I guess my question boils down to whether we can take advantage of the special form of $\boldsymbol{\underline{J}_f}$ when computing the updates $\boldsymbol{\Delta x}$ in each Newton step, $\boldsymbol{\underline{J}_f}(\boldsymbol{x}_k) \boldsymbol{\Delta x} = - \boldsymbol{f}(\boldsymbol{x}_k)$?