问:
Efficient out-of-place arbitrary rank GPU transpose
Summary: Is there an efficient out-of-place GPU tensor transpose operation that scales as $O(n)$ for tensors with $n$ total elements, regardless of the rank $d$?  The naive algorithm costs $O(dn)$, since it requires $O(d)$ worth of index manipulation per entry.

Here is the current implementation of high rank tensor transpose in TensorFlow:

https://github.com/tensorflow/tensorflow/blob/master/tensorflow/core/kernels/transpose_functor_gpu.cu.cc#L27

It does a flat 1-D loop, deconstructs the index using $O(d)$ divisions by strides, and reconstructs the transposed index.  Is there any index trickery / precomputation that would avoid the $O(d)$ cost per entry?


问:
Generalized eigenvalue with null space
Define $S\in\mathbb{R}^{n\times n}$ as
$$S:=H+Q^\top V^{-1} Q.$$
$H,V$ are positive semidefinite.  Here, $H$, $Q$, and $V$ are large, dense matrices but they are structured:  I can write code for matrix-vector products (e.g. $x\mapsto Hx$) but computing $V^{-1}$ or the map $x\mapsto V^{-1}x$ would require iteration.

Given a sparse, symmetric, positive definite matrix $M$, what is the best way to find the $k$ smallest eigenvalues of $S$ satisfying $Sx=\lambda Mx$?

Note that this problem is equivalent to the expanded generalized eigenvalue problem
$$
\begin{pmatrix}
H &Q^\top\\Q  & -V
\end{pmatrix}
\begin{pmatrix}
x \\ z
\end{pmatrix}
=
\lambda
\begin{pmatrix}
M & 0 \\ 0 & 0
\end{pmatrix}
\begin{pmatrix}
x \\ z
\end{pmatrix},
$$
where $z$ is an auxiliary variable.  This is nice because it's symmetric and does not contain $V^{-1}$, so I was hoping to use machinery for generalized eigenvalues.  But the right-hand side of this generalized problem is not invertible, which seems to prevent the use of usual eigenvalue iteration.

ContentLicense: CC BY-SA 3.0

Title: Generalized eigenvalue with null space


问:
Help formulating/finding the general class of this problem
Imagine a bus serving a line with N stations. Each station, $i, i=1,…N$,  has $s_{ij}$ passengers that want to board the bus to go to $j$, $\forall j \neq i$. (one direction). So there are $\sum_j s_{ij}$ passengers waiting at station $i$ to board the bus. Now, suppose that station $m$, is a strategically important  station and we want to make sure there will be enough room on the bus to board waiting passengers. Assume the bus can choose how many people to board at each station, regardless of their destination. We now want to determine how many passengers at each of the stations (prior to station m) should board the train to ensure everyone at station $m$ can board the bus.

What type of optimization problem is this? Typical network optimization problems involve maximizing flow or capacity, but not this sort of problem. How would one go about modeling this? Are there any examples of problems similar to this?


问:
Numeric integration over Dirac delta
I'm trying to solve the following integral numerically.
$$H(y) = \int dx \, f(x) \, \delta(g(x,y)).$$
For this I chose a representation of the delta-function and employ
convergence with respect to $\eta$
$$\sum_\limits{x_i} \, f(x_i) \, \frac{\eta}{g(x_i,y)^2+\eta^2}. $$

I calculate this for every $y$ such that I can plot $H(y)$ in the end.
The problem is, that I would expect for $\eta \rightarrow 0$ the best results, but there seems to be a line (approx. $\eta=10^{-4}$) upon which I get a bunch of delta-like peaks in $H(y)$ instead of a smooth curve (for e.g. $\eta=10^{-2}$) which I want and would expect.

Title: Numeric integration over Dirac delta


问:
Closed-form Jacobian of SE(3) element with 6-degrees of freedom
Let $A$ and $B$ be two rigid transformations in 3D space that transform things from global to local coordinates. Let their relative transformation be expressed by $W=AB^{-1}$. $W$ can also be expressed by an element of SE(3): $\omega=\ln(W)$ where $\ln()$ is defined to be the inverse of the exponential map.

I'd like to know the closed form solutions for computing the 6x6 Jacobians

$\dfrac{\partial\omega}{\partial\alpha}$ and $\dfrac{\partial\omega}{\partial\beta}$

where $\alpha$ and $\beta$ are the six degrees of freedom of $A$ and $B$ respectively. In plain English: How does $\omega$ change when I wiggle $A$ and $B$? It might be that there's a very simple solution to this and I'm just not seeing it.

ContentLicense: CC BY-SA 4.0


Title: MATLAB ODE15s: imaginary parts

Body:
I have to solve
$$\left\{\begin{matrix}
\ddot{R}=-\frac{3}{2}\frac{\dot{R}^{2}}{R}+\frac{1}{R\rho}\left [P_{v}(T_{\infty})-P_{0}+\frac{(L\rho_{g})^{2}}{hT_{\infty}}\dot{R}-4\mu \frac{\dot{R}}{R}-\frac{2S}{R}  \right ]\\
h=\frac{\lambda}{2R}\left [ 2+\left ( \frac{-6\rho c_{p}\dot{R}}{\pi h}\right )^{1/3}-\left ( \frac{12\rho c_{p}\dot{R}}{\pi h} \right ) \right   ]\\
\ln\left ( \frac{\rho_{g}}{\rho_{c}}\right )=c_{1}\theta^{1/3}+c_{2}\theta^{2/3}+c_{3}\theta^{4/3}+c_{4}\theta^{3}+c_{5}\theta^{37/6}+c_{6}\theta^{71/6}
\end{matrix}\right.$$
where $\theta=1-\frac{1}{T_{c}}\left ( T_{\infty}+\frac{L\rho_{g}}{h}\dot{R} \right )$ and $R(0)=R_{0}$, $\dot{R}(0)=0$.
If a use a pressure $P_{0}=2100$, the code works; if I set $P_{0}=2000$, MATLAB returns this error: "Warning: Imaginary parts of complex X and/or Y arguments ignored."

This is my code:

```<code>function out = densita_20_latente(t,y,R0)
S = 0.073; %N/m
rho = 998;
mi = 1.005e-3;
P0 = 2100; %<-----------------------------------------
L = 2454000;
Tinf = 293;
lambda = 0.5982;
cp = 4184;
Pv_Tinf = 2329.6;
c1 = -2.03150240;
c2 = -2.6830294;
c3 = -5.38626492;
c4 = -17.2991605;
c5 = -44.7586581;
c6 = -63.9201063;
Tc = 647.096;
rhoc = 322;
%
out=[y(2)
    -1.5 * y(2)^2 / y(1) + 1 / rho / y(1) * (Pv_Tinf - P0 +...
    (L * y(4))^2 / Tinf * y(2) / y(3) -...
    4 * mi * y(2) / y(1) - 2 * S / y(1))
    -y(3) + lambda * 0.5 / y(1) * (2 + (-6 * rho * cp * y(2) / pi / y(3))^(1/3)...
    - 12 * rho * cp * y(2) / pi / y(3))
    -y(4) + rhoc * exp(c1*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(1/3)...
    + c2*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(2/3)...
    + c3*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(4/3)...
    + c4*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(3)...
    + c5*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(37/6)...
    + c6*(1-1/Tc*(Tinf+y(4)*L/y(3)*y(2)))^(71/6)) ];
</code>```

run file

```<code>R0=510e-6;
y0=[510e-6; 0; 1172.94; 0.01716];
M = [1 0 0 0; 0 1 0 0; 0 0 0 0; 0 0 0 0];
options = odeset('Mass',M);
[t,y] = ode15s(@(t,y) densita_20_latente(t,y,R0),[0 46.75e-6],y0,options);
t = t*1000000;
y = y*1000000;
[t,y(:,1)];
plot(t,y(:,1))
</code>```

Edit: I am going to simplify the problem I have to solve:

$$\left\{\begin{matrix}
\ddot{R}=-\frac{3}{2}\frac{\dot{R}^{2}}{R}+\frac{1}{R\rho}\left [P_{v}(T_{\infty})-P_{0}+\frac{(L\rho_{g})^{2}}{hT_{\infty}}\dot{R}-4\mu \frac{\dot{R}}{R}-\frac{2S}{R}  \right ]\\
h=\frac{\lambda}{2R}\left [ 2+\left ( \frac{-6\rho c_{p}\dot{R}}{\pi h}\right )^{1/3}-\left ( \frac{12\rho c_{p}\dot{R}}{\pi h} \right ) \right   ]\\
\end{matrix}\right.$$

```<code>function out = convezione_20(t,y,R0)
S = 0.073;
rho = 998;
mi = 1.005e-3;
P0 = 2000; %<-------------------------------------------
L = 2454000;
Tinf = 293;
rho_g = 1/57.79;

lambda = 0.5982;
cp = 4184;
Pv_Tinf = 2329.6;
%
out=[y(2)
    -1.5 * y(2)^2 / y(1) + 1 / rho / y(1) * (Pv_Tinf - P0 +...
    (L * rho_g)^2 / Tinf * y(2) / y(3) -...
    4 * mi * y(2) / y(1) - 2 * S / y(1))
    - y(3) + lambda * 0.5 / y(1) * (2 + (-6 * rho * cp * y(2) / pi / y(3))^(1/3)...
    - 12 * rho * cp * y(2) / pi / y(3)) ];
</code>```

run file

```<code>R0 = 510e-6;
y0 = [R0; 0; 1172.94];
M = [1 0 0; 0 1 0; 0 0 0];
options = odeset('Mass',M);
[t,y] = ode15s(@(t,y) convezione_20(t,y,R0),[0 40e-6],y0,options);
t = t*1000000;
y = y*1000000;
[t,y(:,1)];
plot(t,y(:,1))
</code>```

Same problem as before. With 2100 the code works; with 2000 MATLAB reads "Warning: Imaginary parts of complex X and/or Y arguments ignored." I tried using as exponent $1$ instead of $1/3$ in the expression of $h$, but the warning now becomes "Warning: Failure at t=2.874099e-007.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (1.021085e-021) at time t."


问:
Mass conservation in atmospheric continuity equation numerical solution
My phd project is heavily related to numerical modeling of planetary atmospheres. In particular now I am dealing with a particular expression of the continuity equation, involving a thermodynamic flux. Equations read

$$
\frac{\partial n}{\partial t} = -\frac{\partial \Phi}{\partial z}
$$

$$
\Phi = -nD\left(\frac{1}{n}\frac{\partial n}{\partial z} + \frac{1+\alpha}{T}\frac{\partial T}{\partial z} + \frac{1}{h}\right) -nK\left(\frac{1}{n}\frac{\partial n}{\partial z} + \frac{1}{T}\frac{\partial T}{\partial z} + \frac{1}{H}\right)
$$

Where of course $n = n\left(z,t\right)$ refers to the particular chemical species density the equation solves. $h$ is the atmospheric scale height for the particular molecule while H is the average atmospheric scale height corresponding to the average atmospheric molecule. $D$ and $K$ are the molecular and eddy diffusion coefficients respectively. $T$ is the temperature. $\alpha$ is a parameter related to heat diffusion.

I use a standard implicit method to find a system of linear equations so that its solution yields the density at a new time step. I solve a tridagonal matrix in Fortran basically. For 1º order derivatives I take a 1º order central finite differences and 2º order central finite differences for second order derivatives. At the boundaries I use forward and backward approaches when necessary.

As I stated the solution seems correct when studying the density profiles. But there is something very annoying when at the boundaries fixed flux conditions are applied. Let's see by integrating the equation over all the domain

$$
\int_{z_{bot}}^{z_{top}}dz\left(z,t\right) = \rho \left(t\right) = - \int_{z_{bot}}^{z_{top}}dz\frac{\partial \Phi}{\partial z}
$$

$$
\implies \rho \left(t\right) = - \Phi_{top} + \Phi_{bottom}
$$

Being $\rho\left(t\right)$ the magnitude known as column density. After this long introduction here my doubt comes. I fix $\Phi_{bottom} = \Phi_{top} = 0$ so that we should expect $\rho\left(t\right) \equiv constant $. Nevertheless this is not the case at all, mass is not conserved in my numerical system by very big numbers, integrating for each molecule. The processing of results is undertaken by Python scripts using numpy when necessary.

I have been struggling with this for a month and I draw the following few guesses;

Could be that some parameters are not well treated. For instance there is a point where I need to take the numerical derivative of a big set of magnitudes. When studying the outcome it is very very fluctuating for a smooth curve.

Since the system deals with an atmosphere, the background profile used is a hydrostatic kind one, with densities ranging from $n = 10^{20} cm^{-3}$ to $n = 10^{8} cm^{-3}$. Using double precision numbers in fortran could make that during computation numerical fluctuations at the bottom part of the atmosphere may "mask" or screen what happens at the top, yielding a fake production of particles out of the blue.

This method does not conserve mass.

Am I on the correct way to find an answer?

Thank you for your time and attention.


问:
Preconditioning technique for large sparse non-hermitian matrix
I am attempting to solve a computational acoustics problem that involves solving an underlying sparse matrix. The size of the problem varies with grid size (3D) and fill-in's obviously make direct solution impractical. Important features of the matrix are as follows:

It is non-Hermitian and particularly NOT diagonally dominant.

Regardless of the size of the problem, there is always one row that consists of only one element, far away from the diagonal (and a zero on the diagonal). This is part of the problem's closure condition.

I have attempted to solve this problem using iterative techniques with preconditioning to speed up the solution. Iterative solution without preconditioning yields fairly inaccurate results (by comparison with directly computed results for a small grid). ILU preconditioning works for coarse grids (small matrices ~ 42000 x 42000) but fails for anything over 100,000 x 100,000; keeping a low fill-in factor ends in "Factor is exactly singular" and large fill-in factor ends in one of the sub-matrices being singular during pivoting. I also tried enforcing the one zero diagonal element with machine epsilon before attempting ILU but it still fails.

Relevant problem parameters are as follows:

There are 8 PDE's: conservation of mass, conservation of energy, three momentum equations for interior points of computation (full staggered grid with velocity on faces and pressure, temperature & density on cell centers), and three problem specific momentum equations on certain boundaries.

There are suitable closure conditions applied to the problem to obtain a full ranked system, one of which has been described above (zero diagonal row).

I am using a hybrid scheme, with divergences and certain gradients being constructed by using a finite volume approach, and the rest of the operators being constructed using a finite difference approach utilizing polynomial fitting.

If it helps, I am using scipy sparse libraries to carry out my computations. The iterative method I am using is LGMRES, which works without any problems when provided a suitable preconditioner. Please suggest some viable preconditioning techniques for this kind of problem.

Edit: The sparsity pattern of the matrix is as follows:


问:
Convert scipy integration with one step to matlab integration
Scipy integration allows us to do ode integration one adaptive timestep at a time and do something to it. However, matlab ode needs us to specify a timespan , and determine the adaptive timestep within, but cannot do things within.

From scipy

```<code>ODE.set_initial_value(y0, tlist[0])
numoftimes = len(tlist)
for ii in range(1, numoftimes):
    # while loop up to tlist[ii]
    while ODE.t < tlist[ii]:
        t_prev = ODE.t
        y_prev = ODE.y
        # integrate up to tlist[ii], one step at a time.
        ODE.integrate(tlist[ii], step=1)
        #do something to the output
        ...
</code>```

However, in matlab, I failed to do this, mainly step=1 is not available.
Which ODE solver in Matlab allows me to advance in just one timestep only

Attempt is to make the tspan extremely small in matlab ode suite. But this is then not adaptive. In some region the original ode timestep is very big (non-stiff), but at some regions the ode timestep is very small (stiff).

```<code>tlist = [0 tf]
numoftimes = len(tlist)
for ii = 1:numoftimes
     # while loop up to tlist(ii)
     while t_prev < tlist(ii)
            t_prev = T(end,:)
            y_prev = Y(end,:)
            # integrate up to tlist(ii), one step at a time.
            t    = [t_prev,t_prev+0.0000000000001];
            [T, Y] = ode45(dydt, t, y);
            #do something to the output
            ...
     end
end
</code>```

This then is not good as I have, even if the ode45 does integrate with only 1 step
1. force the timestep to be fixed interval
2. time-consuming as in some non-stiff region the "original adaptive" timestep needs not be that small.

How can I solve this with matlab ode suite?


问:
QR via Householder: less computationally complex variants?
I'm a probabilist and need to do a few computations for a rather big linear least squares problem, so I'm trying to optimize the computation as far as is feasible to me.

In computing the QR factorization of the matrix representation $A$ of an overdetermined system $Ax=b$ using the Householder reflections, the commonly stated algorithms suggest that you compute the outer product $uu^\top$ for some $u \in U^\perp$ and a hyperplane $U$, following a geometric motivation.

This yields the Householder reflector $H = I - \frac{2}{\langle u, u\rangle } uu^\top$.

Now, in order to eliminate subcolumns of the columns of $A$, multiply $A$ from left with $H$.

The computation of $uu^\top$ requires $m^2$ multiplcations.

The direct computation of the product $HA$ requires $m^2n$ multiplications.

The product $AH$ is equivalent (or expands) to $A-\frac{2}{\langle u, u\rangle} u(u^\top A)$ which requires $2mn$ multiplications -- by associativity.

I might be missing something rather obvious. Still, is there a way to avoid this computation on every iteration and reduce the complexity from $2mn$ flops?


问:
Non-cubic blocks in Adaptive Mesh Refinement
Most Adaptive Mesh Refinement softwares/logic have cubic blocks. I have explored the BoxLib library and it ONLY supports cubic blocks. To be more precise it supports non-cubic blocks at the coarsest grid level but fails to maintain this non-cubic size when the boxes are refined. I wanted to know if there are block structured AMR packages/softwares which support non-cubic blocks ?


问:
2D Multilateration with constraints as distance information
In my current problem, I am looking for an algorithm to reconstruct the position of multiple points in the 2D euclidean plane with incomplete distance information, $d_{ij}=||x_j-x_i||\in\text{dom}_a$.
E.g. Information is:

$x_i$: $\{d_{ii},d_{ij}, d_{ik}, d_{il}\}\in \text{dom}_a$

$x_j$: $\{d_{ji}, d_{jj}, d_{jk}\} \in \text{dom}_a$, $\{d_{jl}\} \in  \text{dom}_b$

there are more than 2 domains and all domains are non overlapping.

This "binned distance" information is exact and is available of all points.
I like to find the positions or rather areas with high likelihood, where the points can be found.

Currently I think, I have too less information to solve the nonlinear constrained optimization problem. On the other hand, I believe it should be enough to make a Bayesian approach work.

Information, I also have but currently don't make use of:

the binned distance information are from time series and $\max(dx/dt)$ is known. This might be useful for finding the next solution of a given configuration in time.

Any suggestions how to obtain these regions of possible positions are appreciated.

Edit

To answer the questions I added a picture of the forward problem. Black dots represents all positions I like to reconstruct. For two points (green and blue) I exemplary plotted their bins (circle) for which I get the information.

E.g. For the two points close to the blue center, blue would return:
$d_{\text{blue},1}$ and $d_{\text{blue},2}$ are within bin 2 where as $d_{\text{green},1}$ reports within bin 5 and $d_{\text{green},2}$ within bin 6. I do have the radius information for each bin.

All bins are a subset of the real line. Their union represents the real line. Their intersection is empty.


问:
How to implement the timestepping in TR-BDF2 using matlab ode solver as end-points (upper bound)?
Even though i know that i cannot advance the solver ode23tb, which uses TR-BDF2 itself, in just one timestep, i want to set an upper bound on the integration period after which I do something to output and feed the modified output into ode solver again. The upper bounds serve for limiting the timesteps taken in the integration interval tspan = [t_prev ti[j]], (see below), to be as low as possible, hopefully just one. To implement it, i try to do:

```<code>    y_prev = rho0;
    t_prev = 0
    ti = [list of timestep in  TR-BDF2]
    j = 1;
    while t_prev < tf
         #do something to the output
          ...
         # integrate up to adaptive ti in  TR-BDF2 even though inside the interval ode45 has its own timestep
         tspan = [t_prev ti[j]]
         j = j + 1;
         [T, Y] = ode23tb(dydt, tspan, y_prev);
         t_prev = T(end,:)
         y_prev = Y(end,:)
         ....
    end
</code>```

But, how can I know what ti is beforehand? I look at https://math.la.asu.edu/~gardner/TRBDF.pdf
, but do not have an idea on how to do this.


问:
preconditioning LBFGS?
I want to minimize an energy of the form
$$V_1(\mathbf{x}) + V_2(\mathbf{x})$$
where $V_1$ is much stiffer than $V_2$. When I try to use LBFGS, convergence is extremely slow, as the solution oscillates wildly in the "stiff directions" while searching for the optimum in the non-stiff directions.

Is there a way I can modify how I am formulating the problem to avoid this problem? $V_1$ is nonlinear so the stiff directions change during optimization, though from the structure of the problem I can identify them easily at any given iteration.

Title: preconditioning LBFGS?


问:
Determine Lagrange nodal variables of a simplex $T$
Consider a simplex $T$ in $R^d$ with $N_1(T) = \left\{N_i\right\}_{i=0}^{d}\subset P_1^{*}(T)$ be the Lagrange nodal variables (or nodal evaluation). By the Riesz representation theorem, there exist functions $\lambda_{i}^{*}\in P_1(T)$ (which is the representation of $N_i\in P_1(T)$) such that $N_j(\lambda_{i}) = \int_{T} \lambda_{i}\lambda_{j}^{*} = \delta_{ij}$ for each $0\leq i\leq d$. Show that $\lambda_{i}^{*} = \frac{(1+d)^2}{|T|}\lambda_{i} - \frac{1+d}{|T|}\sum_{j\neq i} \lambda_j$

My attempt: I was thinking that $\lambda_i^{*}$ plays a role like a Fourier coefficient, so I multiply both sides by $\lambda_j$ and integrate over $T$. Then LHS would be $\delta_{ji}$, while RHS $= \frac{(1+d)^2}{|T|}\int_{T} \lambda_{i}\lambda_{j} - \frac{1+d}{|T|}\sum_{j\neq i}\int_{T} \lambda_{i}\lambda_{j}$. Now, for $i=j$, then $LHS=1$, and I was thinking of $\int_{T} \lambda_{i}\lambda_{i} = $  the area of region $T$, but then if that's the case, the RHS would be $(1+d)^2,$ which is certainly not equal to $1$ (contradiction). So I don't think my interpretation is correct, but I could not see how to proceed further without knowing the explicit formula for $\lambda_i$.

My question: Could anyone please help me with this problem? It seems to be such a beautiful identity, yet quite hard to prove. Any thoughts would really be appreciated.


问:
computational tool for higher order Lagrangian interpolation for finite element
In finite element, I can calculate the Lagrangian interpolation shape functions for each degree of freedom in an element, from the the number of nodal degrees of freedom and the number of nodes required in the element, if I know the data points ($x$,$y$ coordinates) required to generate linear algebraic conditions, which is in number of total degrees of freedom in the element.

In my problem, for a 2D element I have 6 degree of freedom in each node, namely $\phi, \phi _x, \phi _y, \phi _{xx}, \phi _{yy}$ and $\phi _{xy}$. As the name indicate, the five degrees of freedom are derivative of the $\phi$ degree of freedom.

If I want to calculate, shape functions for degree of freedom in a quadrilateral of four nodes, I would select a incomplete polynomial of degree five and $24$ terms. But the complexity will increase, if the number of nodes are more.

I would like to know that is there any programming library available which can be used for generation of these shape functions for higher nodes?


问:
Simple steady-state advection problem: do I need FVM with upwind scheme?
I have a 2D (x,y) scalar advection problem that describes net blowing snow ($q$) transport at a point. This takes the form

$$q = A - F*\nabla\cdot(q {\bf \hat u}),$$

where A ($kg\cdot m^{-2}\cdot s^{-1}$) is a spatially variable scalar source term, $F$ is fetch (m) and $\hat u$ unit normal parallel to to the wind direction (unitless); $\hat u$ is spatially variable. Initial conditions are [ I believe] $q(x,y)=0$ (no blowing snow, develop steady state) and 0 flux boundary condition. This equation is from the literature, but feels a bit wonky so I'm trying to piece it together.

I would like to solve this via an upwind-scheme -- blowing snow transport moves down wind and should only transport snow to cells downwind.

Further, I would like to solve this on a unstructured triangular mesh.

My hang up (and it's likely a dumb one):

do I discretize this via FVM and then, for each edge, apply an unwind scheme?

Ah the moment I'd like to keep it simple (can move to a better upwind method later), and just find the 'upwind' triangle (via dot-product signs at the edge?) and compute the difference with the current cell. Will that work? If not what is the best way to piece this together?

I appreciate any help, I feel like this is a silly question, but I can't quite wrap my head around the order in which to do this.


问:
Computing Algebraic Riccati inequality
In my Robust model control, I have got a couple of quadratic Riccati inequalities which need to be solved numerically on MATLAB. My question, Is there function on MATLAB can solve quadratic Riccati inequality or I have to transform it to Linear matrix inequality(LMI) using Schur complement approach and then apply it on MATLAB?. Thank you a lot. </p>

Title: Computing Algebraic Riccati inequality


问:
Diverged HDG solution for 2D incompressible Navier-Stokes test case at SMALL time step. Why?
I wrote a hybridizable Discontinuous Galerkin code for transient Navier-Stokes flow,
with thanks to Martin Kronbichler and Scott Miller's step-51 code in Deal.II. The main algorithm
is from Nguyen and Peraire's JCP paper(2012).

$$\frac{1}{\nu}q_{ij}+\frac{\partial u_j}{\partial x_i}=0\;\;in \;\Omega$$
$$\frac{\partial u_j}{\partial t}+\frac{\partial}{\partial x_i}(q_{ij}+p\delta_{ij}+u_i u_j)=f_{j}\;\;in \;\Omega$$
$$\frac{\partial u_j}{\partial x_j}=0\;\;in \;\Omega$$
$$n_i(q_{ij}+p\delta_{ij})|_{\Gamma_N}=g^N_j$$
$$u_j|_{\Gamma_D}=g^D_j$$
$$u_j|_{t=0}=u^0_j$$
The weak form I've taken is as follows,
$$(e_{ij},\frac{1}{\nu}q_{ij})_{T_h}-(\frac{\partial e_{ij}}{\partial x_i},u_j)_{T_h}+<n_ie_{ij},\hat{u}_j>_{\partial T_h}=0$$
$$(v_j,\frac{\partial u_j}{\partial t})_{T_h}+(v_j,\frac{\partial}{\partial x_i}(q_{ij}+p\delta_{ij}))_{T_h}-(\frac{\partial v_j}{\partial x_i},u_i u_j)_{T_h}+<v_j,n_i\hat{u}_i \hat{u}_j>_{\partial T_h}+<v_j,S(u_j-\hat{u}_j)>_{\partial T_h}=(v_j,f_j)_{T_h}$$
$$-(\frac{\partial r}{\partial x_j},u_j)_{T_h}+<r,n_j \hat{u}_j>_{\partial T_h}=0$$
$$(1,p)_{K}=\bar{p}A_K,\;\;\forall K\in T_h$$
and globally,
$$<\mu_j,n_iq_{ij}+pn_j+S(u_j-\hat{u}_j)>_{\partial T_h}=<\mu_j,g^N_j>_{\Gamma_N}$$
$$<1,\hat{u}_jn_j>_{\partial T_h}=0$$
$$(1,\bar{p})_{T_h}=0\;\;(if\; \Gamma_N=\emptyset)$$
Now my code was tested with the 2D lid-driven cavity case, up to Re=1000. All
tests passed with steady mode. But once the time derivatives are turned on,
the code only converged to the steady result with large time step,e.g. dt=100.
If the time step is small, such as 1.e-2, the solution diverged. In all cases,
the stablization coefficient tau is kept constant as 5.0 .

Anyway, the less the time step, the more unstable the code. This really contradicts
againt our common understanding of instability and convergence.

For simplicity, the time marching was done by backward Euler scheme. The Picard
iteration was used for linearization.

I am just wondering if there's any other aspect in the stability of HDG method I should
pay attention to.

Any possible suggestion will be appreciated.


问:
Stability of the explicit MacCormack Scheme to solve the Navier Stokes equations with Wilcox's K-Omega Turbulence Model
I am solving turbulent pipe flow with an explicit MacCormack scheme and Wilcox k-omega model. The laminar version of the code had three distinct stability criteria which worked fine after implementation.

Now, the turbulent version seems to have stability problem. I tried to upgrade the stability criteria for the new set of equations, but it seems that solution eventually blows up. I have noticed the omega values grow very large through all the domain.

Does anybody know if there are new stability conditions with k-omega model, or I should merely extend MacCormack stability to the new equations?

Also I appreciate it if anyone can introduce a paper or book that discusses two-equation models stability.


问:
How to optimally choose points for multivariable Hermite interpolation?
I have a multi-variate, continuous function $f$ from $R^n$ to $R$, which I can query for its output for any input.
I would like to create interpolation polynomial for it.

In one-dimensional case choosing points is relatively straight forward.
I have error estimation:
$$
f(x)-H(x)={\frac  {f^{{(K)}}(c)}{K!}}\prod _{{i}}(x-x_{i})^{{k_{i}}},
$$
and to choose the interpolation points, I need to solve corresponding minmax problem, as I'm looking to minimize the maximum error for a given number of interpolation points. I'm willing to go with any class of  multivariate polynomials. I'm even willing to go with something other that polynomials if helps me interpolate my function, takes into account derivatives and provides theoretical guarantees (for example, in form of error bound).

In multi-variate case I cannot find such error bounds in the literature nor I unable to derive them myself. Any advice on how should I pick the interpolation points for this function?


问:
nonlinear boundary condition
I have a nonlinear system of algebraic equations with a nonlinear boundary condition of the form
\begin{equation}
f'''(x,t) - (f'(x,t))^3 - \alpha(t) f'(x,t) = 0,
\end{equation}
at $x =0$, where $\alpha(t)$ is constant in space but evolves in time. I use sided finite differences to approximate it and Newton's method to solve the nonlinear system approximated using finite differences. The solution does not converge when the mesh size decreases.

Do you suggest any other method to approximate the bc?

Title: nonlinear boundary condition


问:
Number of particles that could be used in a simulation
For a computer with 8 GB of memory, which is the order of magnitude for the number of particles that could be used in a simulation.

Is the answer $N=8/6\times 10^9$? or do the particles take up more memory than 6N(3N position +3N velocity)?


问:
Finite difference method for coupled PDEs: optimizing performance (time step, iterations per step)
I'm solving coupled PDEs using finite difference method: Incompressible Navier-Stokes and the divergence-free induction equation (Maxwell's equations) with non-uniform electrical conductivity. The dimensionless equations are:

$1) \partial_t \mathbf{u} + \nabla \bullet (\mathbf{u}\mathbf{u}) = -\nabla p + Re^{-1}\nabla^2 \mathbf{u} + N \mathbf{j} \times \mathbf{B}, \qquad \mathbf{j} = Re_m^{-1} \nabla \times \mathbf{B}$

$2) \partial_t \mathbf{B} + Re_m^{-1} \nabla \times (\bar{\sigma}^{-1} \nabla \times \mathbf{B}) = \nabla \times (\mathbf{u} \times \mathbf{B})$

$3) \nabla \bullet \mathbf{u} = 0, \qquad \nabla \bullet \mathbf{B} = 0$

Where $Re,N,Re_m$ are the Reynolds, Interaction and magnetic Reynolds numbers. Also, $\bar{\sigma}$ is the dimensionless electrical conductivity.

The induction equation becomes stiff / ill-conditioned as $\bar{\sigma}$ approaches zero. In my case, $\bar{\sigma}$ is very small and, as a result, I resort to treating $Re_m^{-1} \nabla \times (\bar{\sigma}^{-1} \nabla \times \mathbf{B})$ implicitly in time using diagonally Preconditioned Conjugate Gradient Method (PCG). I use the limiting stopping criteria among 3 choices:

$N_{\text{max-iterations}}$

$\frac{||Ax_k-b||}{||Ax_0-b||} \le \varepsilon_{\text{relative}}$

$||Ax_k-b|| \le \varepsilon_{\text{absolute}}$

Here, $A=I + \Delta t Re_m^{-1} \nabla \times (\bar{\sigma}^{-1} \nabla \times)$. $b = \Delta t (\mathbf{B}^n + \nabla \times (\mathbf{u}^n \times \mathbf{B}^n))$. $x_k$ is the $k$th estimated solution of $\mathbf{B}^{n+1}$ using PCG. Finally, $N,\varepsilon$ are the maximum allowable PCG iterations and tolerances (relative to initial guess and absolute respectively).

Experimenting with stopping criteria and different time steps ($\Delta t$), it seems that there is a non-trivial trade-off in picking $\Delta t$ and stopping criteria to optimize performance, which I'll define as

$\mathcal P = \frac{\text{convective time per step}}{\text{wall clock time per step}} = \frac{\text{convective time}}{\text{wall clock time}}$.

Ultimately, I'm trying to reach a steady solution, however, the coupled equations are non-linear, stiff and can (depending on $Re,N,Re_m$) have seemingly long time scales. If I choose $\Delta t$ too large, or perform too few iterations of PCG on Eq. 2), the solution diverges.

Questions:

1) Are there any methods to optimize the choice of stopping criteria or $\Delta t$ in order to optimize $\mathcal P$?

2) Are there any other tricks I can use to reach steady state faster?


问:
Integral over reference element in $1$D FEM: how to map the quadrature points?
The following is related to a question a asked a few days back 1, but now I would like to focus on just one part of the problem.

I have problems computing the integral over the reference element:

$$ \int_{0}^1 (f - \Pi_h^1(f))^2 dx. $$

I think I understood how to use the basis functions on the reference element and what the values of $f$ should be for the interpolant. The problem is with the first $f$ in the integrand. How can I map it? If the affine map is $x = x_j + \gamma (x_{j+1} - x_j)$ I can substitute it into $f$ and obtain

$$ f(x(\gamma)) = cos(2 \pi (x_j + \gamma (x_{j+1} - x_j))). $$

How can I continue? I was following these beautiful notes on page $49$ but that is very simple.

Using Python to do this I ended up with this:

```<code>def func_ref(z, x, a, b):
    cos_a = np.cos(2*np.pi*x[a])
    cos_b = np.cos(2*np.pi*x[b])
    return np.power((cos_a - cos_a * (1 - z) + cos_b * z) * (x[b] - x[a]), 2)
</code>```

where x is simply the vector containing the nodes, a and b would be the left and right neighbor, z the variable of integration. The first term, $\mathbf{cos_a}$ should be the first $f$ of the integrand but I do not know/understand how to do this. Thank you.

Possible solution:

Is it correct to write

$$ \int_{0}^1 (cos(2\pi (x_i + \gamma (x_{i+1} - x_i)) - cos(x_i)(1-\gamma) + cos(x_{j+1})\gamma)^2 (x_{i+1} - x_i) d \gamma $$

and integrate with respect to $\gamma$?


问:
Trying to find a correlation in data from 2 sensors
I've been trying to find a way to correlate data sets from two sensor types. In theory these measure the same signal. One is in a known format expressed in meters while the other is in numeric form read from a 24 bit value. If anyone could help me find a way to find the approximate meter value of these integer values I would be so appreciative for the assistance.

These sensors may have different offsets but the difference between each set of two measurements should be close. I am listing the known values first then the integer representation from the unknown sensor.

Sample 1

0.3940 meters = 38884

0.3860 meters = -348332

Sample 2

0.2125 meters = 32230

0.2055 meters = -97977

Sample 3

0.6510 meters = 21108

0.7495 meters = -363721

I have many more sample sets. Thanks.


问:
Global optimization with known distributions of some variables
I'm solving simple single-objective multidimensional global optimization problem using various stochastic algorithms like Monte-Carlo, GA and other evolutionary approaches. The task is formulated as follows:

search space $X \subseteq R^n$

objective function $f: R^n \rightarrow R$

must find $\textbf{x}^*=\{x_1^*, ..., x_n^* \} \in X: f(\textbf{x}^*)=\min_{\textbf{x} \in X} f(\textbf{x})$

However, I also have some preferences of $x_i \in [a;b]$ which are presented in set of normalized values (from $0$ to $1$). This set can be interpolated (here no matter how) to probability distributions. So, for given $x_i$ exists function $g_i(x)$. Figure below shows typical distribution.

This doesn't mean that $x_i$ always possess preferable values and in peaks exists smaller value of objective. However, better values of $f$ tends to be near peaks.

The question is how to use this information about $x_i$ preferable values during optimization process and acheive better results (e.g. rate of convergence)? Any recommendations/books/articles with similar situations will be appreciated. I even don't know in what math area I should look for.

I have tried to use simulated annealing — just added $g_i$ (multiplied by $-1$) to objective function with a high weight at the beginning and then gradually reduce it to zero. With this simple scheme I didn't acheive much better results.


问:
elliptic equation with exponential coefficient
I'm trying to solve the following equation
$$\dfrac{\partial}{\partial x}\left(e^{au}\dfrac{\partial u}{\partial x}\right) = 0$$

Of course, this equation can be solved analytically. I am trying to understand the way of working with exponential coefficients.

For the approximation, I use the finite volume method which leads  to

$$F_{i+1/2} - F_{i-1/2} = 0,\qquad F_{i-1/2} = \left.e^{au}\dfrac{\partial u}{\partial x}\right|_{i+1/2} $$

For the further linearization, I'm trying the following two methods:

simple iteration, where $e^u$ is calculated based on the previous iteration

Newton's method (with the series expansion, dropping terms that are proportional to $\delta u^2$ and smaller)

$$u^{k+1} = u^k + \delta u$$

The first one worked well until some value of $a$, and the second one didn't work at all (I have not found the mistake yet)

Is there any way to deal with exponential coefficient?  Will be glad for any advice. Thanks!


问:
Software suggestion for simulating the stacking of sedimenting rods
I'd like to simulate the following problem:

Define some square, horizontal and hard area onto which cylinders
(aspect ratio R) sediment.

Preferably, the source of cylinders
sequentially generates N of these cylinders from a large distance
above the plane. There is gravity.

The cylinders are hard, if possible deformable (bending), non-interpenetrable, and friction / damping forces are large to stabilize the system (avoid infinite bouncing).

Essentially I need a physics simulation software which handles collisions of many-body systems. The number of cylinders would be around 100'000.

I am in fine interested in the stacking dynamics of the cylinders. A measure of this would be the z-position of the center of mass of the sedimented mesh of cylinders in function of N/A (aerial density) and the cylinder aspect ratio.

I have knowledge of Python, and can work with scripting Blender, and should be able to find my around making scripts (C,Python,...) work on GPUs.
Blender apparently incorporates Bulletphysics. The Molecular Add-on for Blender that I tried however doesn't have S.I values of physical units (Gravity is mapped between 0 and 1, I don't know what this is), and I couldn't get the simulation to converge.

Any suggestion would be appreciated !


问:
Understanding MATLAB's `fmincg` optimization function
I'm researching numerical optimization. Recently I've come across a variant of a conjugate gradient method named fmincg.

The function is written in MATLAB and is used in the famous Andrew Ng's course on Machine Learning on Coursera. According to the copyright notice, fmincg was written by Carl Edward Rasmussen.

You can find the full text of fmincg on the MATLAB file exchange.

I'm trying to figure out the algorithm. The structure is quite straightforward as it closely follows the reference text I'm currently using (Numerical Optimization, Nocedal & Wright).

However I'm having a harder time figuring out the details of the line search procedure.
Specifically, here's the code for the line search loop:

```<code>while 1
while ((f2 > f1+z1*RHO*d1) | (d2 > -SIG*d1)) & (M > 0)
  limit = z1;                                         % tighten the bracket
  if f2 > f1
    z2 = z3 - (0.5*d3*z3*z3)/(d3*z3+f2-f3);                 % quadratic fit
  else
    A = 6*(f2-f3)/z3+3*(d2+d3);                                 % cubic fit
    B = 3*(f3-f2)-z3*(d3+2*d2);
    z2 = (sqrt(B*B-A*d2*z3*z3)-B)/A;       % numerical error possible - ok!
  end
  if isnan(z2) | isinf(z2)
    z2 = z3/2;                  % if we had a numerical problem then bisect
  end
  z2 = max(min(z2, INT*z3),(1-INT)*z3);  % don't accept too close to limits
  z1 = z1 + z2;                                           % update the step
  X = X + z2*s;
  [f2 df2] = eval(argstr);
  M = M - 1; i = i + (length<0);                           % count epochs?!
  d2 = df2'*s;
  z3 = z3-z2;                    % z3 is now relative to the location of z2
end
if f2 > f1+z1*RHO*d1 | d2 > -SIG*d1
  break;                                                % this is a failure
elseif d2 > SIG*d1
  success = 1; break;                                             % success
elseif M == 0
  break;                                                          % failure
end
A = 6*(f2-f3)/z3+3*(d2+d3);                      % make cubic extrapolation
B = 3*(f3-f2)-z3*(d3+2*d2);
z2 = -d2*z3*z3/(B+sqrt(B*B-A*d2*z3*z3));        % num. error possible - ok!
if ~isreal(z2) | isnan(z2) | isinf(z2) | z2 < 0   % num prob or wrong sign?
  if limit < -0.5                               % if we have no upper limit
    z2 = z1 * (EXT-1);                 % the extrapolate the maximum amount
  else
    z2 = (limit-z1)/2;                                   % otherwise bisect
  end
elseif (limit > -0.5) & (z2+z1 > limit)          % extraplation beyond max?
  z2 = (limit-z1)/2;                                               % bisect
elseif (limit < -0.5) & (z2+z1 > z1*EXT)       % extrapolation beyond limit
  z2 = z1*(EXT-1.0);                           % set to extrapolation limit
elseif z2 < -z3*INT
  z2 = -z3*INT;
elseif (limit > -0.5) & (z2 < (limit-z1)*(1.0-INT))   % too close to limit?
  z2 = (limit-z1)*(1.0-INT);
end
f3 = f2; d3 = d2; z3 = -z2;                  % set point 3 equal to point 2
z1 = z1 + z2; X = X + z2*s;                      % update current estimates
[f2 df2] = eval(argstr);
M = M - 1; i = i + (length<0);                             % count epochs?!
d2 = df2'*s;
end
</code>```

I would like to understand the following:

before the loop, the point 3 is initialized as: f3 = f1; d3 = d1; z3 = -z1; and the comment says Save point 0. But why does he set z3 = -z1?
 (in case someone's wondering, no, it's - apparently - not a mistake because 1. the algo works and 2. I've tried running it with z3 = z1 and it breaks)

how the author is using the Wolfe conditions and its constants?

interpolation / extrapolation. Where do the equations he's using come from? E.g.:

A = 6*(f2-f3)/z3+3*(d2+d3);
 B = 3*(f3-f2)-z3*(d3+2*d2);
 z2 = -d2*z3*z3/(B+sqrt(B*B-A*d2*z3*z3));

what's going on in the if-else checks that involve the limit variable? Note that before entering the while 1, he set limit = -1

EDIT:

added one more question to the bullet list


问:
Which numerical scheme should be used?
Trying to find a way to solve exponentially nonlinear elliptic equation with complex source term i manage to know about such schemes like Godunov, Lax-Friedrich, MUSCL. I still seraching for literature about them.

The main question is could this schemes do something with
$$ \dfrac{\partial}{\partial x}\left(ce^{\varphi_n-\varphi_c}\dfrac{\partial\varphi_n}{\partial x}\right) = R$$

This schemes traditionally for parabolic/hyperbolic equation, but i could add some time derivative (i erase it because i need stationary solution). Or it is possible to use them with elliptic one? (why not)

I also know that it could possible if only the scheme with it initial condition will be stable. For example, MUSCL scheme deal with discontinueties, large gradient etc. (from Wiki)

Is it possible when $ce^{\varphi_n-\varphi_c}$ change very fast (1e18 in a small area)? Could you adviced me some literature?

P.S. sorry for my english, i working on it.

Title: Which numerical scheme should be used?


问:
Relaxing a variable in MIP
I have this MIP optimization problem, with couple of binary variables; however when I relax one of the binary variables the optimal solution of the objective does not change. But the solving time increase drastically. How could this happen? I noticed that this only happens for just when only I relax that specific variable, not for the rest. Does it like this that when you relax a variable the solution space does not change? is the some MIP optimization logic or property behind this?

Thank you for your help.

Title: Relaxing a variable in MIP


问:
Degree of freedom for elastic wave propagation problem
I am solving a elastodynamics (vector valued elastic wave) equation.

I create the 2D mesh in Gmsh discretised into triangular elements of second order. Therefore, it is my understanding that the element is a triangle with total six nodes, i.e., 3 vertices and 1 node on each of the 3 edges. I have the following questions:

Q.1: The dof for a simple scalar problem is equal to the no. of nodes on the mesh. Therefore, the dof for my problem should be equal to twice the no. of nodes on the mesh right? This is because each node can only experience displacement in 'x' or 'y' direction as a result of wave propagation. Is this correct?

Q.2: Each triangular element on the mesh will have 12 dof. Right?

I read somewhere that the velocity of wave in both directions will also constitute 2 dof, but that doesn't seem to make sense.


问:
How one could choose the value of viscous coefficient for obtaining stable solution of Burgers' equation?
Burgers' equation is a fundamental PDE used in various fields such as number theory, gas dynamics, heat conduction, elasticity, etc. It is crucial especially for developing numerical models for turbulence since nonlinear viscous Burgers' equation can be considered as an approach to the Navier-Stokes equations. It can be stated as
$$\frac{\partial u}{\partial t} + \frac{\partial}{\partial x} \left(\frac{u^2}{2}\right) = \nu \frac{\partial^2 u}{\partial x^2}$$ where $a<x0$ with appropriate initial and boundary conditions given. (it can be determined later)

Let assume that the finite difference numerical method are considered to obtain approximate solution of this viscous equation. In my perspective, there are two decisions need to be made to have stable solutions

What kind of finite difference method can be considered for inviscid and viscous part (or both of them) i.e., Lax-Wendroff, Lax Ricthmyer, Crank-Nicholson, Lax-Friedrich's, Roe's scheme ?

How one can decide the value of the artificial diffusion number ($\nu$) ? Should we consider its relation with Reynolds number since it is related with the Navier-Stokes equations ?

Note : This viscous Burgers' equation is not going to be considered as linearized. It will be considered in a nonlinear form.


问:
Fast solution of a heptadiagonal linear system
I have a linear system of the form $\mathbf{A}\mathbf{x} = \mathbf{f}$. If the length of the vector $\mathbf{x}$ is $N$, meaning that there are $N$ unknowns, then the matrix $\mathbf{A}$ has seven bands and the matrix is sparse. The matrix is not symmetric. It looks like this:

$\pmatrix{a_1&a_2&a_3&a_4&&&&&&\\
b_1&b_2&b_3&b_4&b_5&&&&&\\
c_1&c_2&c_3&c_4&c_5&c_6&&&&\\
d_1&d_2&d_3&d_4&d_5&d_6&d_7&&&\\
&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&&\\
&&&&&&\\
&&&k_1&k_2&k_3&k_4&k_5&k_6&k_7&&&\\
&&&&l_1&l_2&l_3&l_4&l_5&l_6&l_7&&&\\
&&&&&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&\ddots&&\\
&&&&&&\\
&&&&&&&w_1&w_2&w_3&w_4&w_5&w_6&w_7\\
&&&&&&&&x_1&x_2&x_3&x_4&x_5&x_6\\
&&&&&&&&&y_1&y_2&y_3&y_4&y_5\\
&&&&&&&&&&z_1&z_2&z_3&z_4\\
}$

and the vector $\mathbf{f}$ looks like this:
$\pmatrix{f_1\\f_2\\f_3\\f_4\\\vdots\\\\f_i\\f_{i+1}\\\vdots\\\\f_{N-3}\\f_{N-2}\\f_{N-1}\\f_{N}}$

I know that the general method is to use a LU decomposition and as the matrix is constant for me, I can pre-compute its decomposition and store it and use it to solve for the unknowns. But I would like to know if there is any method to solve such systems faster than the LU decomposition. Methods may not be restricted to direct methods only.


问:
Is there any implementation of Nesterov's paper: Gradient methods for minimizing composite functions
I worked on a method based on this paper, "Gradient methods for minimizing composite functions". But unfortunately, every time I code it doesn't work. Even I investigate inequalities of the paper numerically, but I couldn't figure out what's the problem. Is there any matlab implementation of the paper available? I search the internet but I didn't find any thing? In my work, I don't need to find Lipschitz constant as it did in paper. </p>


问:
A least square problem with a fixed mean constraint and a subspace constraint
Let $V_1,\ldots,V_n$ be $n$ vector subspaces of a Hilbert space, $y_i\in V_i$ for each $i$ and $\overline{x}$ be a fixed vector. I want to solve the optimization problem:
\begin{equation*}
\begin{aligned}
& \underset{x_1,\ldots,x_n}{\min}
& & \sum_{i=1}^n \| x_i - y_i\|_2^2 \\
& \text{subject to}
& & \frac{1}{n}\sum_{i=1}^n x_i = \overline{x},\\
& && \ x_i \in V_i \text{ for }i=1,...n.
\end{aligned}
\end{equation*}
I don't have a basis for the $V_i$ a priori but I can use the projection operators $P_i$ as I please. The setting is supposed to be an infinite dimensional Hilbert space but any finite dimensional help is appreciated too. The $V_i$ are not necessarily orthogonal nor does their direct sum equal the whole space necessarily.

If the $V_i$ are orthogonal and direct sum into the whole space, I believe the solution is $x_i = y_i + N*P_i(\overline{x} - \overline{y})$.

Feel free to add in any assumptions you might need.

Asked on MSE originally, sorry if reposts are frowned upon.


问:
A stable method for solving monontoe HJB equation
I am considering solving HJB equation of the form
$$
v_t=g(a(x)v_x),\quad x\in \mathbb{R}, t>0,
$$
with initial condition $v(0)=v_0$. Here $g:\mathbb{R}\to \mathbb{R}$ is Lipschitz and monotone increasing. I hope to construct a scheme stable in $L^\infty$ norm.

Here is what I tried. For the linear $g$, i.e. $g(p)=cp$ with possible space dependent $c>0$, I consider the upwind and downwinding scheme；
$$
D_iv^{n}=\begin{cases}\frac{v_{i+1}^{n}-v_{i}^{n}}{h},& a_i>0,\\
\frac{v_{i}^{n}-v_{i-1}^{n}}{h},& a_i<0 \end{cases},
$$
and update the solution by
$$v^{n+1}_i=v^{n}_i+\tau c_ia_iD_iv^n.$$
Then for $a_i>0$, we have
$$v^{n+1}_i=v^{n}_i+\frac{\tau c_ia_i}{h}(v_{i+1}^{n}-v_{i}^{n})=(1-C_i)v_i^n+C_iv^{n}_{i+1}\le |v^n|_\infty,$$
with $C_i=\frac{\tau c_ia_i}{h}>0$, while for $a_i<0$, we have
$$v^{n+1}_i=v^{n}_i+\frac{\tau c_ia_i}{h}(v_{i}^{n}-v_{i-1}^{n})=(1-C_i)v_i^n+C_iv^{n}_{i+1}\le |v^n|_\infty,$$
with $C_i=-\frac{\tau c_ia_i}{h}>0$.

However, for Lipschitz and monotone increasing $g$, the scheme becomes
$$v^{n+1}_i=v^{n}_i+\tau g(a_iD_iv^n).$$
To perform similar analysis, I write
$$g(a_iD_iv^n)=g(a_iD_iv^n)-g(0)+g(0),$$ thus
$$|v_i^{n+1}|=|v_i^n+\tau (g(a_iD_iv^n)-g(0))|+\tau|g(0)|,$$
then using Lipschitz $g$ I have
$$|v_i^{n+1}|=|v_i^n|+L\tau |a_iD_iv^n|+\tau|g(0)|.$$

But the similar analysis can't process since I don't know how to write $|v^{n}_i|+C_i|v^{n}_{i+1}-v^{n}_{i}|$ as a convex combination. Also I haven't used the monotone increasing assumption.

I have also tried the scheme discussed here, but I still can't demonstrate the stability due to a similar reason.

Is my simple upwinding scheme a stable scheme to solve this problem? I am not sure whether the proof fails due to my poor techniques or the scheme is not stable in general for this type of $g$. Is there any other scheme commonly used for this problem?


问:
Variable elimination in linear programming
I have a linear program of the form

$$\underset{P,\;g}{\text{Minimize}}\hspace{3mm}c^Tg$$
\begin{align}
\hspace{17mm}\text{Subject to}\hspace{3mm}AP_{\cdot,j}&=\begin{bmatrix}
-g\\
d
\end{bmatrix}\\
P_{ii}&=0\\
|P_{\cdot,j}|&\leq P_{\max}\\
g&\leq G_{\max}\\
g&\geq 0
\end{align}

with variables: $g\in{\bf R}^k$, $P\in{\bf R}^{n\times n}$,

and problem data: $c\in{\bf R}^k$, $A\in{\bf R}^{m\times n}$, $d\in{\bf R}^{m-k}$, $P_{\max}\in{\bf R}^n$, $G_{\max}\in{\bf R}^k$.

I'm only interested in $g$, and thus I was wondering if there was a way to somehow eliminate the variable $P$, thereby significantly reducing the dimension of the search space.

If $A$ is invertible it's almost certainly doable.

If $A$ has a left-inverse it might be doable.

If $A$ has more columns than rows and full row-rank it might be doable.

If $A$ has more columns than rows without full row-rank then I don't think it's doable, however this is the case I'm most interested in.


问:
Solving an integral equation in Python
I have to solve the following equation for $x(i), 0 \leq i \leq 1$:

$$ y(i) = x(i)^{-a} \int_0^1 y(j)x(j) dj \left(\int_0^1 \mathcal A(j) x(j)^{1-a} dj\right)^{-1} \int_i^1 \left( \int_0^x x(j)^{1-a} dj \right)^\frac{a+b}{a-1} A_i^{b - 1} dx $$

where $A_i$ and $y(i)$ are known and I use the short cut

$$ \mathcal A(j) = \int_j^1 A_z^{1 - b} dz$$.

I guess there is no point in trying to solve it symbolically with Sympy. How would I go on solving this numerically using Python? Do I discretize $x(i)$ into 100 or 1000 points and hope to brute force find a solution to the system of equations? Or is there a more clever approach?

Title: Solving an integral equation in Python


问:
How to implement boundary conditions on Finite Difference WENO5 scheme for the Euler equations
I'm implementing a Finite Difference WENO5 with Lax-Friedrich flux splitting on a uniform, structured grid to solve the 2D Euler equations of fluid dynamics on a rectangular domain in cartesian coordinates.

$$\frac{\partial U}{\partial t}+\frac{\partial F}{\partial x}+\frac{\partial G}{\partial y}=0$$

where:

$$U = (\rho,\rho{u},\rho{v},\rho E)^t$$
$$F = (\rho{u},\rho{u^2}+p,\rho u v, \rho uH)^t$$
$$G = (\rho{v},\rho u v,\rho{v^2}+p, \rho vH)^t$$

$$E = \frac{p }{\rho (\gamma-1)} + \frac{1}{2}(u^2+v^2)$$
$$H = \frac{\gamma p }{\rho (\gamma-1)} + \frac{1}{2}(u^2+v^2)$$

Although the interior scheme seems to work fine, I am having trouble with the numeric implementation of boundary conditions.

The WENO5 interpolation is done on the fluxes, to approximate the derivative on the node $i$ as: $$\left( \frac{\partial F}{\partial x}\right) _{i}=\frac{f^{WENO}_{i+1/2}-f^{WENO}_{i-1/2}}{\Delta x}$$ and then update the conservative variables' vector via a Runge-Kutta method.

I want to use ghost nodes, and I know there has to be 3 ghost nodes on each end of the domain (because of the stencil of the WENO5) but I do not know what values I have to assign to the ghost nodes, or what kind of extrapolation I have to use to find those values. In my case the physical boundaries coincide with the leftmost and rightmost (lowermost and uppermost for the other direction) points of the grid, respectively.

These are the three types of boundary conditions I'm working with at the moment, namely:

Reflecting wall (x-direction): $$\begin{bmatrix}\frac{\partial \rho }{\partial x}=0 \\ u=0\\ \frac{\partial v }{\partial x}=0 \\ \frac{\partial p }{\partial x}=0 \end{bmatrix} $$

Reflecting wall (y-direction): $$\begin{bmatrix}\frac{\partial \rho }{\partial y}=0 \\ \frac{\partial u }{\partial y}=0 \\ v=0\\ \frac{\partial p }{\partial y}=0 \end{bmatrix} $$

Dirichlet condition with fixed value on all primitive variables: $$\begin{bmatrix} \rho =\rho_0 \\ u =u_0 \\ v =v_0 \\ p =p_0 \end{bmatrix} $$

Thanks in advance for your answers.


问:
Numerically compute PDF given a function
Consider $[0,1]$ with the Lebesgue measure $m$ and $f:[0,1]\to \mathbb{R}$, and $x$ a uniformly distributed random variable in $[0,1]$. Then, $f(x)$ itself define a new random variable.

We can then define the Cumulative density function of $f$ as $$P_f (y) :\,= m(f^{-1} (-\infty, y)) \, ,$$ and its PDF as $$\rho _f (y) :\,= \frac{dP_f (y)}{dy} \, , $$ for all $y\in {\rm range} (f)$.

A straightforward numerical computation of the probability density function (PDF) of $f$ involves inverse derivatives, or, conversely, $1/f'$. As numerical differentiation is numerically ill-conditioned, when $f'$ becomes very low, this is numerically problematic. I'm looking for one of two things:

An out of the box function in Matlab that computes a PDF for a given function $f$.

A numerical Algorithm for this task.

Thanks


问:
Order of local and global truncation error in a finite volume scheme
Can the order of the local truncation error be higher than the order of global truncation error in a finite volume method?</p>


问:
Proportional equality constraints
Consider a node $s$. Let's assume that there are three outgoing arcs from node $s$ namely $(s,i)$,$(s,j)$ and $(s,k)$. Corresponding to each of these arcs, there is a flow proportion value $t_{sj}\in (0,1]$. It means that the amount of flow along the arc $(s,j)$ should be $t_{sj} *f_{ps}$ in which $f_{ps}$ is flow over the single arc going to node $s$. Now in my LP I have the following constraints regarding node $s$:

$f_{ps}  =  f_{si}+f_{sj}+f_{sk}$

$ f_{si} =  t_{si} *f_{ps}$

$ f_{sj} =  t_{sj} *f_{ps}$

$ f_{sk} =  t_{sk} *f_{ps}$

In a special case, just assume that $\sum_j {t_{sj}}=1$. For example assume that $t_{sj}$ are probability values.

Now, my question is that how I can getting rid of these proportional equality constraints, i.e. constraints  2-4 in our given example. I am looking for a way to modify the structure of the network in such a way to include these constraints implicitly.

Title: Proportional equality constraints


问:
Connectivity and Clustering using Eigenvectors and the Fiedler Vector
Going off of the answer here:
sorting adjacency matrix by the Fiedler vector

So here, Jesse the answerer plotted the first 3 eigenvectors associated with non-zero eigenvalues of the Laplacian against each other in order to see more distinguishable clusters.

My 2 questions about this are:

1) Why did he choose specifically the first 3 eigenvectors to plot, as opposed to for example, just the first 2?  Like what determines how many eigenvectors you should compare?

2) Is there a way to quantify connectivity between clusters?  For example, when looking at the graph Jesse originally made, is there a way to tell from plotting the eigenvector values that the top cluster is more closely connected to the second top most cluster than the bottom most cluster?


问:
Optimization based integration for MPM
I'm considering implementing (just for simplicity) the unconstrained implicit optimization based integration for Material Point Method as described in Chenfanfu Jiang's thesis on MPM (the minimization algorithm starts on page 101 – "A.3.1 Unconstrained minimization").

I have difficulties understanding how would be the algorithm implemented.

Are the statements below (in)correct?

The goal of the algorithm is obtaining new grid velocities in response to external and internal forces of the material represented by material points.

The function to be minimized is $E(\mathbf{v}) = \sum_i\frac{1}{2}m_{\mathbf{i}} \lVert \mathbf v_i - \mathbf {v_i}^n \rVert + \Phi(\mathbf {x_i}^n + \Delta t\mathbf v_i)$. (from Jiang's thesis)

It's derivative is $\mathbf h(\mathbf v) = \mathbf M\mathbf v - \Delta t \mathbf f(\mathbf x^n + \Delta t \mathbf v) - \mathbf M\mathbf v^n$.

Based on the thesis I assume the approach to minimizing $E$ is by the Newton method. The Newton method requires inverse of hessian $\mathbf{\mathit{H}}$. However we're in fact interested in $\Delta \mathbf v$ that would minimize (the linearization of) $E$, so we can instead solve $\mathbf{\mathit{H}}\Delta \mathbf v=\mathbf h$.

When solving the system via Conjugate gradient, we only need to know the product $\mathbf{\mathit{H}}\Delta \mathbf v$, which can be thought of as differential $\delta \mathbf h = \mathbf M\Delta \mathbf v - \Delta t^2 \delta \mathbf f(\mathbf x^n + \Delta t \mathbf v) = \mathbf{\mathit{H}}\Delta \mathbf v$.

From [Stomakhin et al. MPM snow paper] we know how to compute $\delta\mathbf{f_i}$ (page 6 – "6 Stress-based forces and linearization") – i.e. force differential for the grid velocity at node $\mathbf i$.

My assumptions:

When solving $\mathbf{\mathit{H}}\Delta \mathbf v=\mathbf h$ via Conjugate Gradient (step 4. in the algorithm)

${\Delta \mathbf v}$ is a column vector consisting of component-wise grid velocity differentials – $[{\Delta v}_{1x}, {\Delta v}_{1y}, {\Delta  v}_{1z},\,\,\, {\Delta v}_{2x}, {\Delta  v}_{2y}, {\Delta v}_{2z}, \,\,\,\ldots]$ (length of $\Delta \mathbf v$ is therefore $3n$, where $n$ is the number of grid nodes (containing velocity))

When multiplying $\mathbf{\mathit{H}}\Delta \mathbf v$, we're in fact always processing 3 rows at the same time - i.e. we process one grid node at time (which consists of x,y,z velocity components).

When computing residual $\mathbf r = \mathbf h - \mathbf{\mathit{H}}\Delta \mathbf v$, we can compute its consecutive 3 components corresponding to grid node $\mathbf i$ (like in assumption 1.1.) by writing $\mathbf{r_i} = \left(m_{\mathbf i}\mathbf v_{\mathbf i} - \Delta t \mathbf f_{\mathbf i}(\mathbf x^n_{\mathbf i} + \Delta t \mathbf v_{\mathbf i}) - m_{\mathbf i}\mathbf {v_{\mathbf i}}^n\right) - \left(m_{\mathbf i}\Delta\mathbf v_{\mathbf i} - \Delta t^2 \delta \mathbf {f_i}(\mathbf x^n_{\mathbf i} + \Delta t \mathbf v_{\mathbf i})\right)$.

When computing $\mathbf r \cdot \mathbf r$, I ignore the fact that each 3 consecutive elements of vector $\mathbf r$ correspond to 1 grid node, and just compute the dot product like with any other vector.

Additional question:

Is the line search (step 7 of the algorithm) necessary?

Title: Optimization based integration for MPM


问:
sum of n numbers with the power k
I am searching for a good algorithm to solve problem 487 on project-euler. I dont want code or something like that, I only want the name of the algorithm thats best suited. Till now I think of either using Faulhabers Formula or using the stirling numbers of second kind to compute the sum.
https://projecteuler.net/problem=487
If there is a better algorithm please let me know.</p>

Title: sum of n numbers with the power k


问:
How does this Constrained Minimization algorithm work?
I don't fully understand the subsection 3.2 Constrained minimization of this paper. In particular, I don't understand the first step "Register active set" and the definition of projection $P(x)$.

Simplified algorithm (stripped-down version of the paper's one):

Select initial guess and project it: $x^{(0)} = P(x'^{(0)})$

Loop until precise enough

Register active set

Solve $\nabla^2E \cdot \Delta v = \nabla E$ via Conjugate Gradient to obtain $\Delta v$

Use line search to find step $\alpha$

Take step: $v^{(n+1)} = v^{(n)} + \alpha\Delta v$

Project $v^{(n+1)}$

Projections. Let $P(x)$ be the projection that applies $P_{bp}$ to
  $x_p$ for all body-particle pairs $(b, p)$ that are labeled as active
  or are violated ($\phi_b(x_p) < 0$).
  Note that pairs such that
  $\phi_b(x_p) = 0$ (as would be the case once projected) are considered
  to be touching but not violated. The iterates $x^{(i)}$
  obtained at the end of each Newton step, as well as the
  initial guess, are projected with $P$.

...
  We begin by representing our collision objects (indexed
  with $b$) by a level set, which we denote $\phi_b$ to avoid confusion
  with potential energy. By convention, $\phi_b(x) < 0$ for points $x$
  in the interior of the collision object $b$. Our collision constraint
  is simply that $\phi_b\left(x_p^{(n+1)}\right) \gt 0$ for each simulation
  particle $p$ and every constraint collision object $b$. With such a
  formulation, we can project a particle at $x_p$ to the closest
  point $x'_p$ on the constraint manifold using $x'_p = P_{bp}(x_p) = x_p - \phi_b(x_p)\cdot\nabla\phi_b(x_p)$. ...

Register active set. Let $E'$ be the objective that would be computed in the unconstrained case. The objective function for constrained optimization is $E(x) = E'(P(x))$.

  Compute the gradient $\nabla E'$. Constraints that are touching ($\phi(x_p) = 0$) and for which $\nabla E'\cdot\nabla\phi \ge 0$ are labeled as active for the remainder of the Newton step. All others are labeled as inactive. No constraint should be violated at this stage. Note that $E'(x^{(i)}) = E(x^{(i)})$ is true before and after every Newton
  step, since constraints are never violated there.

The emphasized sections are the ones I don't understand.

From the above I assume that the projection function $P(x)$ is defined as
  $$\begin{align}
\phi_b \left( x_i \right) &= \text{the shortest distance from point $x_i$ to the surface/boundary of $b$ } \\
P \left(x_i\right) &=
\begin{cases}
  x_i - \phi_b \left( x_i \right) \cdot \nabla\phi_b \left( x_i \right)
    & \,\,\, \text{if $\phi_b \left( x_i \right) < 0$ or when particle $i$ is active} \\
   x_i
    & \,\,\, \text{otherwise}\\
\end{cases}
\end{align}
  $$
  As can be seen from the algorithm, and from the fact $E(x) = E'(P(x))$, the projection function has tangible effect only on such particles that are "active" (since no constraint is violated inside of the loop).

What is the purpose of labeling a pair $(b, p)$ as "active"? And why are labeled only those pairs for which $\nabla E'\cdot\nabla\phi \ge 0$?

I do currently understand this as selecting those particles $p$ that violated the constraints and were thus projected on the level-set manifold (or just happened to be at the precise boundary of the manifold) and at the same time the direction of force (represented by $-\nabla E'$) is trying to push the particle outside the manifold ($\nabla\phi$ is the surface normal of the level-set pointing inside the manifold at the position of particle $p$).

Do I understand it correctly that the projection function has effect only on the derivatives of the objective function $E'(x^{(i)}) = E(x^{(i)})$ (because otherwise it degrades to identity)?

Why exactly is $P = x_i - \phi_b \left( x_i \right) \cdot \nabla\phi_b \left( x_i \right)$ for the "active" particles? By definition, for the "active" particles $\phi_b \left( x_i \right) = 0$, so the value of $P$ is $x_i$ – the same as for the inactive particles.


问:
Numerically inverting an exponentially growing function (defined by Chebyshev polynomials)
Assume a function $M(t)$ strictly increasing, essentially growing exponentially, and asymptoptically  growing at a known rate $\bar{g}$, i.e. $\lim_{t\to\infty}M'(t)/M(t) = \bar{g}$

In a set of awful integro-differential equations, I need to use the inverse of $M(t)$, denoted by $q(\cdot)$ for simplicity, i.e. $t = M^{-1}(M) \equiv q(M)$.

The function $M(t)$ will be solved numerically with some sort of spectral collocation method with a polynomial basis (see below).  My general problem is: with a guess on $M(t)$, how can I find $q(M)$?

Since this equation is growing exponentially, I don't think I can use a polynomial basis for $M(t)$ directly.  Instead, I suspect that I should guess a Chebyshev basis for $g(t)$ constraining for some large $T$ that $g(T) = \bar{g}$.  If so, then
$$M(t) \equiv \exp\left(\int_0^t g(\tau) d \tau \right)$$

Which works fine for my computations.  In particular, it is easy with a Chebyshev basis to find the partial integrals $\int_0^t g(\tau) d \tau$ for $t$ nodes of the Chebyshev basis (through Calculating integrals for a function approximated by Chebyshev polynomials)

But if so, how can I find the $q(M)$ function, or a basis of it in $M$ space?  Obviously, I could just evaluate $M(t)$ at the Chebyshev nodes of $t$ to find the appropriate $M$, and then use some sort of interpolation to find the approximating $q(M)$, but I would prefer an approach that has a chance of using auto-differentiation by throwing it into my collocation method.

For example, is there a way to translate this into a nonlinear ODE using the inverse-function theorem, and if so does this simplify substantially given that this is a Chebyshev basis?


问:
Finding a permutation that makes a matrix lower triangular
I have a system of linear equations in form of $AX=b$ where $A_{m\times n}$, $X_{n\times 1}$ and $b_{m\times 1}$. Coefficient matrix $A$ is quite sparse. However, using a practical LP solver like LINGO it is clear that after permutation of rows, it turns out to be like a lower triangular matrix. I do not know what the permutation function should be. Since the matrix dimension is not square, I cannot use the LU decomposition to solve the system efficiently. Can you please let me know an efficient method for linear systems with non-square coefficient matrix?</p>


问:
Comparing the solutions to a multi-objective optimization problem
Suppose I have a multi-objective optimization problem, and I wish to find solutions using two different methods/algorithms. The result of each algorithms is a Pareto front. Comparing two different methods for a single-objective optimization problem is simple: just see which one produces the better result for the single objective. However, for a multi-objective problem, one would need to compare the resulting Pareto fronts, supposedly.

My question then is: is it possible to compare two Pareto fronts? How would you say that one algorithm performs "better" than the other?

Thanks, and any insight is valued.


问:
Rhie and Chow Pressure Velocity Coupling
In a collocated grid, which velocity is used in the convective term in the momentum equation? Is the Rhie and Chow constructed face velocity or an average of the adjacent cell center values(in a structured grid) used?</p>


问:
Shallow water equations with moving body
I was wondering if there are any shallow water equation solvers out there that can include rigid body motion for the surface the water flows on?</p>


问:
How to test for convergence (smoothness of Pareto front) in DEAP
In the DEAP algorithms (see documentation here), I notice that we need to specify the the number of generations (NGEN). I was advised that convergence has been achieved if the Pareto curve is smooth. How can I monitor for convergence in DEAP?

Update:
It is possible to monitor for convergence by specifying a "smoothness" value in the stats. However, I'm still confused about how to define "smoothness". For example, consider the Knapsack problem specified here. How do we monitor for smoothness in this example?


问:
Closed form PDF/CDF using Orthogonal Polynomial Expansion (gPC)
Consider a random variable which is given by an orthogonal polynomial expansion in one parameter (or polynomial chaos expansion PCE), i.e. ,
$$ f(\alpha) = \sum\limits_{n=0}^{\infty} \hat{f} (n) \psi _n (\alpha) \, , $$

where $\alpha$ is uniformly distributed in $[-1,1]$,  $\psi _n (\alpha)$ are the orthogonal Legendre polynomials $$\hat{f} = \frac{1}{2}(f,\psi _n )_{L^2[-1,1]}  \, .$$

My question: is there a closed form formula or an approximation for the cumulative distribution function (CDF) or the probability distribution function (PDF) of $f$ in a closed form, using only $\psi _n$ and $\hat{f} (n)$?

This is somewhat related to this post, and was asked by me before on MO here


问:
Good numerical method for solving the Kadomtsev Petviashvili equations. Is there an analytical solution?
I need to solve the  Kadomtsev Petviashvili (KP) equations

$$\partial_x(\partial_t u+u \partial_x u+\epsilon^2\partial_{xxx}u)+\lambda\partial_{yy}u=0 $$

where $$\lambda=\pm 1 \;.$$
My questions are:

1) Which is the best numerical method to solve them? Or maybe a good method that is simple?

2) What does wikipidia mean with: "the KP equation is completely integrable"? I don't think it means KP has a general analytical solution, otherwise why are people still coming out with numerical methods for solving it? See for example this paper.


问:
First approximation to the TSP in a non-complete Graph
I'm trying to solve the Travelling Salesman Problem in a non-complete graph $(G,E)$ using genetic algorithms.

My problem is that I can't find a good first approximation by the usual greedy algorithms, as long as I can't assure an arbitrary edge $(u,v)$ to be in $E$. I tried to add to $E$ new edges (with almost infinity weight) in order to make it complete, bit then it turns to be way more inefficient.

I can't find a good solution on internet. So can anyone help me?


问:
Stochastic gradient descent for large deterministic optimization problems
The Wikipedia page for SGD describes optimizing a function $f = \sum f_i(\theta;x_i)$ by successively approximating gradients from random subsets of the data, while most literature poses the problem as optimizing $f = \text{E}[F(\theta,\zeta)]$ for some random variable $\zeta$. I can see how these are asymptotically equivalent when $\zeta$ is a sample from all available observations and $F(\theta; \zeta) = \sum_{x_i\in \zeta}f_i(\theta;x_i)$, and intuitively, it seems that the later formulation is more descriptive of the actual optimization task. However, it isn't apparent to me that these problems necessarily have the same solution. Under what conditions are the solutions equal (convexity, etc.)? Are there any theoretical guarantees that a stationary point of the stochastic approximation be in some $\epsilon$-neighborhood of a stationary point of the deterministic function when these assumptions are not met?

EDIT: Wikipedia says the method converges to a global optima for convex and pseudo convex functions. I'm curious if there is any established error analysis between a deterministic solution and it's stochastic approximation, or are these algorithms used because they make solutions for normally intractable problems possible, and so error analysis is generally impossible/unimportant? I'm trying to decide if an SGD, or one of several online L-BFGS methods, is a flexible enough framework for a tool where I expect a large range of dataset sizes (from 10 to 10s of thousands) and dimensions, but is still generally too many dimensions for most second order methods, without losing accuracy or significant speed when the problem is small enough for other deterministic methods.


问:
Numerical Double integration with endpoint singularity in scipy Python gives incorrect answer
I am trying to integrate the following function in Python,

$\int_{0}^{\infty}\int_{0}^{\infty} \dfrac{e^{-x-y}}{B(x,y)}dx dy$,

where $B(x,y)$ is the beta function - $B(x,y) = \int_0^{1}a^{x-1}(1-a)^{y-1}da$.

As the above function was not directly integrable in scipy.integrate, I did a change of variable, $x\rightarrow e^{-x},y\rightarrow e^{-y}$ to have

$\int_{0}^{1}\int_{0}^{1} \dfrac{1}{B(-\log (x),-\log (y))}dx dy$.

scipy.integrate.dblquad(lambda y,x: 1/ scipy.special.beta(-math.log(x),-math.log(y) ),  0, 1, lambda x: 0, lambda x: 1)

returns $6.81\ldots$ with error in order of $10^{-6}$. But  integrating the above functions (with and without change of variable) in Mathematica (using NIntegrate) returns $2.0505\ldots$. I do not understand this huge difference in answer when python returns an answer with $10^{-6}$ accuracy. Whats happening? Is there a better way to integrate it in python? I plan to perform the integration many number of times as part of an algorithm written in python.

$\frac{1}{B(-\log (x),-\log (y))}$ looks like this

Edit

Can someone with Windows 10 and Anaconda reproduce the results I get?

I get the incorrect answer using Windows 10, Python 3.5.2 |Anaconda 4.3.0 (64-bit), scipy version 0.18.1. I am not sure about the BLAS linkage.

I get the correct answer in Linux using Python 2.7.12, scipy version 0.18.1


问:
linear relaxation of an optimization problem
I'm reading an article lately, and there is one point which confuses me. So, we have the following constrained binary quadratic problem.

min $x^{T}Qx$

with the constraints that $Ax≤b$

and $x\in {0,1}^{n}$

This can be relaxed into the following:

min $\sum_{1\leq i <j \leq n}2q_{ij}y_{ij} + \sum_{i=1}^{n}q_{ii}x_i$

with the constraints that

$Ax\leq b$

$y_{ij}\geq x_i+x_j-1 \text{ for any } i,j \text{ with } i<j \text{ and } q_{ij}>0 $

$y_{ij}\leq x_i \text{ for any } i,j \text{ with } i<j \text{ and } q_{ij}<0$

$y_{ij} \leq x_j \text{ for any } i,j \text{ with } i<j \text{ and } q_{ij}<0$

$0\leq x_i \leq 1\text{ for any } i=1,2,...n$

$y\geq 0$

So, the solution of the relaxed optimization problem would serve as a lower bound of the optimal solution for the original optimization problem. Could anyone help me explain the relaxation part? I was thinking to prove

$\sum_{1\leq i <j\leq n}2q_{ij}y_{ij} + \sum_{i=1}^{n}q_{ii}x_i \leq \sum_{j=1}^{n}\sum_{i=1}^{n}x_iq_{ij}x_j$ where $x_i$ is either 0 or 1.

It seems that the term $\sum_{i=1}^{n}q_{ii}x_i$ of the left part would cancel out the term $\sum_{i=1}^{n}x_i^{2}q_{ii}$ of the right part.
Then, what shall I do next? Any comments would be greatly appreciated.


问:
Unclear definition of objective function for MPM
While reading this paper on MPM (page 48), the definition of the objective function $E(\boldsymbol v_{\boldsymbol i})$ didn't make sense to me, although I completely understand the definition of (what I suppose to be) it's derivative $\boldsymbol h(\boldsymbol v^{n+1})$.

Excerpt from the paper:

We will show that solving

$$\boldsymbol h(\boldsymbol v^{n+1}) = \boldsymbol M \boldsymbol v^{n+1}  -\Delta t \boldsymbol f(\boldsymbol x^n + \Delta\boldsymbol v^{n+1}) - \boldsymbol M \boldsymbol v^n = \boldsymbol 0$$

is equivalent to minimizing the following objective function:

$$E(\boldsymbol v_{\boldsymbol i}) = \sum_{\boldsymbol i} \frac{1}{2} m_{\boldsymbol i}\lVert \boldsymbol v_{\boldsymbol i} - \boldsymbol v_{\boldsymbol i}^n \rVert^2 + e(\boldsymbol x_{\boldsymbol i}^n + \Delta t \boldsymbol v_{\boldsymbol i}).$$

I believe $\boldsymbol v^{n+1}$ and $\boldsymbol v^{n}$ are velocities of all grid nodes of cartesian grid, whereas $\boldsymbol v^{n+1}_{\boldsymbol i}$ and $\boldsymbol v^{n}_{\boldsymbol i}$ are the velocities of one particular grid node.

The questions:

Why is $E(\boldsymbol v_{\boldsymbol i})$ taking only single velocity vector $\boldsymbol v_{\boldsymbol i}$ as its parameter ($\boldsymbol v_{\boldsymbol i}$ = the velocity of the grid node at position $\boldsymbol i$) while it's summing over all grid nodes $\boldsymbol i$ (the sum $\sum_{\boldsymbol i}$) – isn't that nonsense? Shouldn't it take all grid nodes' velocities $\boldsymbol v$ instead?

Isn't it nonsense to sum $e(\boldsymbol{\hat x})$ $i$-times, given that:

$e(\boldsymbol{\hat x}) = \sum_p V^0_p \Psi(\boldsymbol{\hat F}_p(\boldsymbol{\hat x}) )$ - i.e. sum over all particles $p$ (page 45/equation 191)

supposing $E(\boldsymbol v_{\boldsymbol i}) =   \sum_{\boldsymbol i} \left(\frac{1}{2} m_{\boldsymbol i}\lVert \boldsymbol v_{\boldsymbol i} - \boldsymbol v_{\boldsymbol i}^n \rVert^2 + e(\boldsymbol x_{\boldsymbol i}^n + \Delta t \boldsymbol v_{\boldsymbol i})\right)$ (note the parentheses)

I'd suggest to define $E$ instead as
$$E(\boldsymbol v) = e(\boldsymbol x^n + \Delta t \boldsymbol v) +  \sum_{\boldsymbol i} \frac{1}{2} m_{\boldsymbol i}\lVert \boldsymbol v_{\boldsymbol i} - \boldsymbol v_{\boldsymbol i}^n \rVert^2$$
Is this redefinition of $E$ correct? If not, why is the original formulation correct?


问:
Optimizing estimator of composed functions when function is known
Note: This was cross-posted from comp-sci, as I didn't know this community existed!

I have a problem which I'm looking to see if there is literature on:

Consider three types of actors, a Director, Aggregator, and Follower. The Director talks to multiple Aggregators, the Aggregator talks to multiple Followers. The director gives simple commands to the Aggregator, who then can portion out the commands to each Follower as she sees fit.

The question then is, how can the Director know whether or not an Aggregator can complete a given request at time step T? T+1? If T is now, then this is relatively easy as the Aggregator can simply sum up the capabilities of each of the Followers and return it, saying "this is what I can do in total." The problem then is, if the Director is doing her own optimization problem in deciding how to direct two different Aggregators, given each command she must be able to predict how the Aggregator will react. I.e., what are the Aggregator's capabilities in time step T+1 given a direction D in time step T?

Due to the complexity of the Aggregator, I think it's likely best to create an estimator across the current state of the Aggregator (or all of their followers). Further, the Aggregator can give whatever knowledge it wants to the Director, though ideally this would be of constant size no matter the size of the set of Followers. This is therefor an optimization problem on the Aggregator's side for creating an optimal estimator, which is well studied.

What I want to know: Is there any literature about creating an estimator when you know the exact function you are estimating? So if I am the aggregator, how can I tell the Director, in a constant space for the amount of information I give no matter the number of followers, the information require to best predict I would react to a given direction?

Put another way: How can I create an optimal estimator in which I have full knowledge of the shape and composition of the functions I am estimating?

I can think of a number of places where this might occur, potentially in networking, data compression, etc., though I think I don't have the necessary vocabulary. I've tried searching through aggregation, optimization, estimation, forecasting, prediction, which a number of different combinations of keywords relating to these. There is a ton of generalized work on making estimators and predictors, though none which target the assumption that they already have perfect knowledge of that which they're estimating.

[Post note after finding this stack exchange]
I read a bit, finding this for instance, and felt that my question was still unique. In this case all functions are known beforehand.


问:
constrained quadratic binary problems and quantum adiabatic evolution
I'm going through an article with title "Solving constrained quadratic binary problems via quantum adiabatic evolution" (reference 1). And there are several points confusing me a lot.

This article is aimed to solve the CBQP (constrained binary quadratic programming) with the following format.

\begin{align}
&\min &x^{T}Qx \cr
&\text{subject to} &Ax\leq b
\end{align}
and $x\in \lbrace0,1\rbrace^{n}$, where $Q\in Z^{n\times n}$ and $A\in Z^{m\times n}$. Let's call this optimization problem the problem $P$.

The outline is like this. Suppose there is a UBQP (Unconstrained binary quadratic programming) oracle, and with the successive application of LP(linear programming), the lagrangian dual of $P$ (or lower bound of $P$ can be provided) can be solved. And then with the branch-bound-approach, the problem $P$ can be solved.

I can understand almost of it until the section 5 counting the solution density on page 9. I'm not exactly sure how the branch-bound process is combined with the optimization process to finally tackle this problem. Could anyone point a direction or share some thoughts? Any comments would be greatly appreciated.

References

Ronagh, P., Woods, B., & Iranmanesh, E. (2015). Solving constrained quadratic binary problems via quantum adiabatic evolution. arXiv preprint:1509.05001.


问:
How to model pedestrian flow through subway systems?
I'm a New Yorker and take the subways every day. I have a growing interest in understanding the distribution of paths people take on the subways to work every day.

I.e. if there are $n$ subway stations, I want a temporal distribution over the paths $\omega = (n_i, n_j, t)$ where $n_i$ and $n_j$ are nodes representing subway stations and $t$ is a continuous variable representing time.

I am given the expected arrival and departure times for each subway station as well as the number of people who enter and exit each station at discrete moments in time.

I am hoping someone can give me ideas on how to translate my deterministic model into one I can fit statistically, as well as reason about theoretically

Ideal Theoretical Model

Let $N$ be the total population in the subways at any given moment in time and let $N_i$ be the number of individuals in station $i$, and $N_ij$ be the number of individuals on the subway line heading from station $i$ to station $j$. The total number of individuals in transit
$$T = \sum{N_{ij}}$$

and the total number of individuals in standby is
$$S = \sum{N_i}$$

Then $$N = T+S$$

I have 1.5 undergrad physics classes under my belt so I've been thinking about the system with the following analogy. Imagine the subway station to be a system of $n$ point masses(each point mass is a subway station) suspended in space, where $N_i$ represents the height of each point mass.

Each adjacent subway station exerts a downwards and upwards force on its neighbors. In particular $F_{ij}$ is the force exerted by $i$ on $j$ such that $$F_{ij} = g(N_{ij})$$ where $g$ is some unknown function.

In this light, the net force on station $i$ is defined as

$$F_{\text{net}, i} =  \sum{ F_{ji} } - \sum{ F_{ij} }$$

In other words, at each moment in time, every station exerts a force on its neighbors. Actually, to be more precise, this force based model implies there are actually $n^2+n$ particles between stations and $n$ stations, therefore a total of $n^2 + 2n$ particles. The velocity of a particle is the rate at which its height decreases and the mass of a particle is a coefficient $m_{ij}$.

Having interpreted the subway station as a group of point particles with pairwise interactions under an unknown force mechanic, we can also define a potential energy configuration.

$U(N)$ defines the potential energy configuration of the system and $U_{ij}(N)$ defines the potential energy of a particle $ij$, such that

$$U(T(t)) + K(S(t)) = E(t)$$

$$K_{ij} = \frac{1}{2}m_{ij}v_{ij}^2$$

$$U_{ij} = -\int^{N_{ij}}_{0}{F_{ij}} = -N_{ij}m_{ij}$$

Fluctuations in $T$ and $S$ represent changes in the internal potential and kinetic energies of the system. Fluctuations in $N$ represents net work being done on the system, i.e.

$$W = \Delta E$$

Issues in the Model

I still haven't defined $g$, nor defined the potential and kinetic energy functions.

Statistical Modelling and my Question

At this point, I have a cute little physical interpretation of the subway system. But I am having trouble translating my dynamical model into one I can statistically fit to find a distribution over the ultimate paths taken by individuals over the course of a day.

In practice, I don't know what "modelling the subway dynamics" will actually mean?

Ought I to fit coefficients to $g$ , $U$, and $K$?

How to formulate a rigorous statistical question given this abstract physical model.

Or

Alternative ways of finding the distribtion of trips, including links to literature that has already solved this problem.


问:
Computing only the $k$ biggest eigenvalues and eigenvectors with Scalapack
Given that there are eigensolvers in Scalapack that use a divide and conquer method, is there any way we can use Scalapack functions to only compute the first $k$ dominant eigenvalues and corresponding eigenvectors?

This is for a spectral clustering code.


问:
Optimisation of matrix exponential
I have a 7000x7000 sparse matrix (scipy), which I want to exponentiate. I've tried using scipy.sparse.linalg.expm, which works quite well for smaller matrices (takes a few seconds for a 1000x1000 matrix) but it takes too long time to compute for the matrix in question. Is there a way to optimise this or work around the time complexity of this algorithm (apparently O(n^2))?
For background, I've tried the same with R package {expm} with similar results, so it's apparently more about the algorithmic approach to the problem rather than the specific software used. Do you have any ideas?</p>

Title: Optimisation of matrix exponential


问:
The right algorithm to predict from selected data samples
I'm trying to come up with the right algorithm for a system in which the user enters a few symptoms and the system has to predict or determine the likelihood that a few selected symptoms are associated with those existing in the system. Then after associating them, the result or output should be a specific disease for the symptoms.

The system is comprised of a series of diseases with each assigned to specific symptoms, which also exist in the system.

Let's assume that the user entered the following input:

```<code>A, B, C, and D
</code>```

The first thing the system should do is check and associate each symptom (in this case represented by alphabetical letters) individually against a data-table of symptoms that already exist. And in cases where the input doesn't exist, the system should report or send feedback about it.

And also, let's say that A and B was in the data-table, so we are 100% sure that they're valid or exist and the system is able to give out the disease based on the input. Then let's say that the input now is C and D where C doesn't exist in the data-table, but there is a possibility that D exists.

We don't give D a score of 100%, but maybe something lower (let's say 90%). Then C just doesn't exist at all in the data-table. So, C gets a score of 0%.

Therefore, the system should have some kind of association and prediction techniques or rules to output the result by judging the user's input.

Summary of generating the output:

```<code>If A and B were entered and exist, then output = 100%
If D was entered and existed but C was not, then output = 90%
If all entered don't exist, then output = 0%
</code>```

What algorithm would be used to produce or come up with the above system concept?


问:
Order of accuracy of FVM discretization
I've recently got interested in CFD and started a small project by solving the radial Reynoldsequation. Why the Reynoldsequation? I recently encountered it through my studies and somehow got stuck :)

I would greatly appreciate it, if someone were so kind to have a read through my steps and point out my thought mistakes or any other errors.

In a first approach, I would like to solve it without the source term, which reduces it basically to a two dimensional diffusive process.

$\frac{\partial}{\partial r}(rh^3\frac{\partial p}{\partial r}) + \frac{1}{r}\frac{\partial}{\partial \theta}(h^3\frac{\partial p}{\partial \theta}) = 0$

If one were to imaging a ring, with two Diriclet Boundary Conditions:
$p(r_I) = 2$ and $p(r_O) = 250$, so I don't need to worry about the singularity at $r=0$. I also have a periodic condition along the angle at $\theta = 0$ and $\theta = 2\pi$
The actual pressure values don't matter, as long as they are larger as zero. Negative pressures don't really work :D In the end, I shall have a unsymmetrical problem with several boundary conditions on the ring, that's why I don't wish to reduce this to a one dimensional one.

Following the approach by Patanker or any other book on Finite Volume Methods (such as The Finite Volume Method in Computational Fluid Dynamics by F. Moukalled), the semi discretization with the midpoint rule will give me a second order accurate approximation:

$r_eh_e^3\Delta \theta(\frac{\partial p}{\partial r})_e - r_wh_w^3\Delta \theta(\frac{\partial p}{\partial r})_w + \frac{\Delta r}{r_i}(h_s^3(\frac{\partial p}{\partial \theta})_s - h_n^3(\frac{\partial p}{\partial \theta})_n)$

That can be confirmed using the Taylor analysis for example.

For convenience, I wrote the code in Matlab, so the East direction is $j+1$ and the South direction is $i+1$ for a matrix given by $M(i,j)$ Lowerbound letters, such as e and s denote the faces midway between two points.

I use a structured grid and equidistant in r and $\theta$ ($\Delta r \neq \Delta \theta$)

Now I'm getting to the part, where I've been scratching my head for some time now.

If one were to approximate the gradient at the faces of the control volume like:
$(\frac{\partial p}{\partial r})_e = \frac{p_E - p_i}{\Delta r}$ and $(\frac{\partial p}{\partial r})_w = \frac{p_i - p_W}{\Delta r}$ and similar for the south and north faces, the formal order of accuracy is what?

Taylor analysis yields a formal order of accuracy of $1$ for the gradient approximation as above, while my books and other sources (internet etc.) refer to the combined discretization (see below) as of order $2$ on a uniform grid and others as of order $1$. They contradict themselves, which is kind of annoying.

$\frac{r_eh_e^3(p_E-p_i)\Delta \theta}{\Delta r} - \frac{r_wh_w^3(p_i - p_W)\Delta \theta}{\Delta r}$

I understand, that when I choose $d_1 = r_eh_e^3\Delta \theta$ and $d_2 = r_wh_w^3\Delta \theta$ (similar for South and North) and in the special case of $d_1 = d_2$, the above discretization would become a second order central difference scheme look a like. But I don't have that case.

$\frac{d (p_E - 2p_i + p_W)}{\Delta r}$

I implemented the discretization and performed an order of accuracy study yielding an order of accuracy of $1$ by comparing the solution of the scheme in radial coordinates to the exact solution. I did that by integrating the resulting pressure field to obtain the hydrostatic lift force. I can now compare the numerical one to the exact one.

Is this actually allowed? The order of accuracy I get is basically 1 (+/- 0.003) for fine grids and starts to oscillate for coarser grids as expected. Would it be better to directly compare the resultant pressure fields?

That's basically it. I sifted through so many texts that I'm now confused to a point, where every explanation I try to come up with is clouded in self doubt.

I would like to find out, which formal order of accuracy I should obtain and  how I can adapt the discretization to achieve an order of accuracy of $2$.

The above discretization's of the gradients are basically a linear profile between two points, but I'm not sure how a higher order profile would look like.

Edit 1:

I use a staggered grid, where the pressure is stored at the centre of the cell and the height at the cell faces. The radii, for example $r_e$, are not approximated but direclty computed at the position of the faces by $(j+c \pm 0.5)*\Delta r$, where c is a constant to account for the inner radius, as my domain does not extend to $r=0$ but starts at $r_I$.

To be able to use the exact solution, the height is currently uniform and could be taken out of the equation, but will be later evaluted through the interpolation of a matrix. Additionally, the density and viscosity could also be added, currently assumed constant, and would also be evaluated at the cell faces.

I also updated the above equations, such that capital letters denote a cell centre and the small letters the cell faces, where i is also at the centre of a cell and conincides with the centre of my stencil.

Edit 2:

I'm currently using these three books:

1: Patanker: Numerical Heat Transfer and Fluid Flow Link

2: F. Moukalled: The Finite Volume Method in Computational Fluid Dynamics Link

3: J.H. Ferziger: Computational Methods for Fluid Dynamics Link

Patanker doesn't go much into detail regarding the order of accuracy, while in the second book they note (page: 216), that the semi discretization is of order two, but the linear profile used to approximate the gradient is of first order, which would make the complete approximation be of order $1$ from my point of view.

The third book goes more into detail regarding the linear profile used for the gradient on pages: 72-74. They note, that the linear gradient approximation is of formal order $1$, but behaves as a second order approximation on a uniform grid, when one chooses the faces midway between two neighbouring points, as I have. On Page 75 it states, that the scheme converges asymptotically in a second order manner.

Several lecture notes on the internet just mentioned, it would be second order
accurate without going much into detail.

\enclose{horizontalstrike}{After going through the texts again, from my understanding, I should obtain a order of accuracy of one, which I do. Yet, I would like to be certain of it^^}

Answer
May I answer my own question? Well, anyway.

I've now been able to obtain an order of accuracy of $2$. I had the error norm computed wrong or in other words, I did not scale the error of the obtained lift force to the grid.

Title: Order of accuracy of FVM discretization


问:
Can convolution be generalized to 2D from several 1D convolutions?
I have a function $\mu(x)$ that I need to convolve with several functions $\nu_n(x)$. I'm currently taking the DFT of $\mu$ and multiplying it individually with the N DFTs of the functions $\nu_n$ and then taking the inverse DFT to obtain the convolution.

Can I instead define
$m(x, y) = \mu(x), n(x, y) = \nu_y(x)$ for discrete $y$, take the 2D DFT to obtain $M(u, v)$ and $N(u, v)$, which I then multiply element wise as $C(u, v) = M(u, v) N(u, v)$ and then take the inverse 2D DFT to obtain the same result $c(x, y)$ as I did elementwise in the above method?


问:
Generate discrete set of points in a feasible region
I have two vectors which specify the bounds $x_{min}$ and $x_{max}$ of the sample space. Also, it has to satisfy the linear constraint $Ax \leq b$.

How to generate an evenly spaced set of points, with some control over grid spacing, which satisfy these constraints and this includes both the boundaries and the interior? In other words, how do I generate a 'mesh' for a feasibility region? The feasibility problem in this case

\begin{equation*}
\begin{aligned}
& \underset{x}{\text{minimize}}
& & 1 \\
& \text{subject to}
& & Ax \leq b \\
&&& x_{min} \leq x \leq x_{max}
\end{aligned}
\end{equation*}


问:
Literature on numerical solving based on multiple meshes?
Consider solving a differential equation system for 1D, 2D, or 3D.

It involves various input and output "field" variables, which, more often than not, correspond to various physical quantities relevant to the problem.

The typical discretization step is to agree on a mesh, and then solve the system. We might refine the mesh in some areas of the spatial-domain (1D, 2D, or 3D) to increase resolution.

However, each "field" has its own high-resolution requiring regions, which are, in general, different from every other field. E.g., considering 1D only, a sin(x) function is relatively flat around pi/2, but steep around 0, whereas cos(x) is flat around 0, but step around pi/2. That means a non-uniform mesh optimized for a sin 1D field is not optimal for a cos field, and vice versa. This means we need a mesh for sin, and another mesh for cos, if we are to do things in the most optimized way.

Despite spending many years in the area of scientific computing, I'm surprised that I've never come across such a discussion that highlights the need for "multi-mesh numerics". Granted it has its own problems, e.g., the proliferation of meshes for intermediate results, computational cost of conversion, and so on, but academically, it's a direction worth pursuing.

So my question is, can anyone help me point to a body of literature that covers this topic? Or is there a different term used for this that I'm not familiar with?