问:
Problem in using N-BFGS-B in optim
I am a beginner user of R. I am trying to maximize log likelihood function with the bounded parameters. The function is a kind of gamma mixture model which try to capture unobserved heterogeneity across individual. I chose the optim's L-BFGS-B given that the function contains two parameters (i.e., alpha and r) which should have positive values. However, the optimization does not take place and returns the initial parameters with the  message below.

```<code>#Sample data
data<-matrix(c(1,36,547,2,54,464,3,92,415,4,114,1106,5,10,1038),ncol=3,byrow=TRUE)
colnames(data)<-c("ID","freq","time")
</code>```

\begin{align}
&\min_{r, \alpha}\; r \log \alpha + \sum_{j=0}^{n_i-1}\log(r+j) - (n_i + r)\log(t_{is} + t_{ic} + \alpha)\\
&\text{such that}\\
&\qquad r > 0\quad \text{(shape parameter)}\\
&\qquad \alpha > 0\quad \text{(scale parameter)}
\end{align}

In the code below I use following terms for function.

ni = data[i,2] ## ni is the frequency of individual i

tis+tic = data[i,3] ## tis+tic is the total time of individual i

```<code>#Likelihood function I want to optimize
ll <- function(theta){
alpha<-theta[1]  #scale parameter >0
r<-theta[2] #shape parameter >0

ll_i=c()
s=c()

for(i in 1:5){

for(j in 0:data[i,2]-1){
s[j+1]<-log(r+j+1)
}

ll_i[i]<-r*log(alpha)+sum(s)-(data[i,2]+r)*log(alpha+data[i,3])
s=c()
}

print(ll_i)
ll<-sum(ll_i)
return(-ll)
}

# Initial parameters
init<-c(alpha=1,r=5)
# L-BFGS-B optimization
ans<-optim(par=init,fn=ll,method="L-BFGS-B",lower=c(1e-10,1e-10),upper=c(Inf,Inf))

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"
</code>```

Any hint towards a successful optimization?

Title: Problem in using N-BFGS-B in optim


问:
Convergence of KKT equations for discrete parameter estimation problems
Consider a discrete constrained optimization problem:

$$
\mathbf{q}_*^h= \arg \min {\cal J}^h(\mathbf{x}^h[\mathbf{q}^h],\mathbf{q}^h)
$$
subject to the (weak-form) constraint
$$
F^h[\mathbf{x}^h(\mathbf{q}^h),\mathbf{q}^h](\mathbf{w}^h) = 0 \; \forall \, \mathbf{w}^h
$$
with $\mathbf{x}^h \in \Omega^h$ (the discrete domain) and $\mathbf{w}^h$ the discrete test functions. We can assume $F^h$ is a consistent and stable discretization of a corresponding continuous problem on the discrete domain $\Omega^h$. The first order necessary conditions that define our local minimum $\mathbf{x}^*$, with $\lambda_*^h$ the discrete adjoint variable (or Lagrange multiplier) are:

$$
F^h[\mathbf{x}^h,\mathbf{q}_*^h](\mathbf{w}^h) = 0 \;\; \forall \, \mathbf{w}^h \\
F_{\mathbf{x}}^h[\mathbf{x}^h,\mathbf{q}_*^h](\mathbf{w}^h,\lambda_*^h) = {\cal J}_\mathbf{x}[\mathbf{x}^h,\mathbf{q}_*^h] \;\; \forall \, \mathbf{w}^h \\
F_{\mathbf{q}}^h[\mathbf{x}^h,\mathbf{q}_*^h](\mathbf{w}^h,\lambda_*^h) = {\cal J}_\mathbf{q}[\mathbf{x}^h,\mathbf{q}_*^h] \;\; \forall \, \mathbf{w}^h
$$

How sensitive is the discrete optimum $\mathbf{q}^h_*$ to consistent perturbations to the primal problem? By consistent I mean a discrete perturbation that vanishes when evaluated at the exact (continuous) KKT triplet $(\mathbf{u},\lambda,\mathbf{q})$:

$$
F^h[\mathbf{x}^h(\mathbf{q}^h),\mathbf{q}^h](\mathbf{w}^h) + \Delta^h [\mathbf{x}^h,\mathbf{q}^h](\mathbf{w}^h) = 0 \; \forall \, \mathbf{w}^h
$$

What bounds do we need to prove on $\Delta^h$ (or its derivative) to guarantee the perturbed KKT system still converges to $(\mathbf x_*^h,\mathbf q_*^h)$? I'm aware the answer may well be problem dependent, so I'd welcome some references to existing literature for discrete model problems (linear diffusion, etc.).

ContentLicense: CC BY-SA 3.0


问:
Reflecting boundary condition posed as a Riemann problem
I am trying to implement a solver for the Euler/Navier Stokes equations. I have a problem implementing boundary conditions for the wall. I am using an unstructured solver.

A lot of literature says that at the wall you simply reflect the normal velocity, and then pose it as a Riemann problem. So let me take the momentum equation. The U here is $\rho u$. The F (flux) is $\rho u^2+p$.
$u$ is velocity, $\rho$ is density, $p$ is pressure. Assume that the grid is such that $u$ is normal velocity.

Now, the way literature guides toward implementing this boundary condition is to take ghost values such that boundary conditions at the wall are satisfied. So reflect the velocity, and copy density, pressure.

The velocity will come to zero at the wall using this. But if you look at the flux, it will not. I have compared Lax-Friedrichs, HLL and I don't see how F becomes $p$ as it should.

Any help?

For reference I will put up this link

http://arc.aiaa.org/doi/abs/10.2514/6.2014-2923


问:
How to obtain values in physical space for a given spectrum?
My question falls under purview of turbulent flows. I want to add an initial perturbation, for which I have a given energy spectrum (say$ E(k)=ak^4e^{-bk^2} $). The steps involved in getting these perturbation is as follows:

Choose random values for $u\in (-1,1)$

Fourier transform:  $u(x,y,z)\rightarrow\hat{u}(k)$

Make $\hat{u}$ solenoidal $\hat{u}'=\hat{u}-\hat{k}.\hat{u}/k^2$

Scale $\hat{u}'$ with the given energy spectrum

Take inverse transform:  $\hat{u}'(k)\rightarrow u(x,y,z)$

How to proceed when I have to give perturbations at discrete points, in which case $u(x,y,z)$ is available at discrete points $(x_i,y_j,z_k)$?


问:
Heat equation from implicit scheme with Neumann B.C
To solve the heat equation from implicit scheme subject to Neumann boundary condition we can write:
$$
T_i^{j+1}-T_i^{j}=\alpha (T_{i+1}^{j+1}-2T_{i}^{j+1}+T_{i-1}^{j+1})
$$
$$
\textbf{A} T^{n+1} = T^n
$$
$$
-\alpha T_{i-1}^{n+1}+(1+2 \alpha)T_i^{n+1}-\alpha T_{i+1}^{n+1} =T_i^n
$$

Initial and boundary conditions are :
$$
T(L/4,t)=100 \quad ,   \frac{\partial T(0,t)}{\partial x} =\frac{\partial T(L,t)}{\partial x} =0
$$

where alpha is $kdt/dx^2$, $k$ is the coefficient of thermal conductivity and L is the length of the rod. from the boundary condition we have :

```<code>T(2,j)= T(0,j)
T(n+1,j)=T(n-1,j)
</code>```

so we have a tridiagonal matrix with main diagonal of ($1+2\alpha$) and subdiagonal and superdiagonal of $-\alpha$. The first and last row of the matrix has changed by Neumann.B.C to

```<code>A(1,1)    A(1,2)  :    1+2\apha    -2\alpha
A(n,n-1)  A(n,n)  :   -2\alpha     1+2\alph
</code>```

Here is my implementation of the problem in FORTRAN:

```<code>program heat_rod
!
! To solve the heat eauation for a rod with
! Tt(x,t) =sigma* Txx(x,t)
! T(x,t) = 0 and T(L/4,t)= 100,
! N.B.C : Tx(0,t) = Tx(L,t) = 0
! x  over [0,1] and t over [0,0.1]
! Variables :
! dl, dd, du subdiagonal, maindiagonal and superdiagonal elements respectively
! To plot the result:
! echo 'splot "data" nonuniform matrix with lines' | gnuplot --persist

implicit none

integer,parameter :: n = 11, m = 31      ! n=21, m=81
integer :: i,j,index_x
real(8) :: alph,l,time
real(8) :: dl(n),dd(n),du(n),d_x,d_t,sigma
real(8) :: T_new(n),T_old(n)

open(unit=7, file="T.txt")

l = 1.d0
time = 1.d-1
d_x = l/real(n-1)
d_t  = time/real(m-1)
sigma = 1.d0 ! coefficient of termal conductivity
alph = (sigma * d_t)/(d_x*d_x)
!print *, alph
if (alph > 0.5) then
print*, "Solution may be unstable"
!stop
endif

du = -alph
dl = -alph

!du(1) = 0.d0         !Dirichlet.B.C
du(1) = -2.d0 * alph  !Neumann B.C
!dl(n) = 0.d0         !Dirichlet B.C
dl(n) = -2.d0 * alph  !Neumann B.C

dd = 1.d0 + 2.d0 * alph
!dd(1) = 1.d0
!dd(n) = 1.d0
!dd(1) = 1.d0 + 2.d0 * alph
!dd(n) = 1.d0 + 2.d0 * alph

T_old = 0.d0
index_x = l/4 *n
!print *, index_x
T_old(index_x) = 100.
!print *,T_old

! print the x positions at the first row
write(7,120) n,((i-1)*d_x,i=1,n)
120 format (I3,999f12.2)

write(7,110) 1,(T_old(i),i=1,n)

do j = 1,m
   call tridag(dl,dd,du,T_old,T_new,n)
   write(7,110) (j)*d_t ,(T_new(i),i=1,n)
   T_old = T_new
enddo

! Printing the results
!write(*,120) ((i-1)*h,i=2,nx-1)
!120 format (7x,' t',7x,'x = ',f5.2,9f11.2)
!do j = 1,m
!    write(7,110) (T_plot(i,j),i=1,n)
!enddo
110 format (f6.3,999f12.6)

!close(7)

 contains

SUBROUTINE tridag(a,b,c,r,u,n)
INTEGER n,NMAX
REAL(8) :: a(n),b(n),c(n),r(n),u(n)
PARAMETER (NMAX=500)
! Solves for a vector u(1:n) of length n the tridiagonal linear set
! a(1:n) , b(1:n) , c(1:n) , and r(1:n) are input vectors and are not modified.
! Parameter: NMAX is the maximum expected value of n .
INTEGER j
REAL(8) :: bet,gam(NMAX)      !One vector of workspace, gam is needed.
if(b(1).eq.0.)pause 'tridag: rewrite equations'
!If this happens then you should rewrite your equations as a set of order N − 1, with u2
!trivially eliminated.
bet=b(1)
u(1)=r(1)/bet
do  j=2,n
!Decomposition and forward substitution.
gam(j)=c(j-1)/bet
bet=b(j)-a(j)*gam(j)
if(bet.eq.0.)pause 'tridag failed'
!Algorithm fails; see below.
u(j)=(r(j)-a(j)*u(j-1))/bet
enddo
do  j=n-1,1,-1
!Backsubstitution.
    u(j)=u(j)-gam(j+1)*u(j+1)
enddo
return
END SUBROUTINE
end program
</code>```

I doubt about the the way I defined the boundary conditions.
How can I validate the simulation results? Thanks in advance for your guide or comment.


问:
Jump conditions on interior boundaries
I can not undrestand when we have jump on interior boundaries, if the solution be continuous. In  the interior boundaries, derivatives have only different unit outward normal vector and this means we can discard all interior derivatives terms. But in some text i see that they are different and the notation of $[\frac {\partial u}{\partial n}]=\frac {\partial u^+}{\partial n}-\frac{\partial u^-}{\partial n}$ is used. Is it related to the test space?
In fact, I am working on adaptive finite element methods for Poisson's equation, based on the book: The finite element method,Theory,implementations and applications  written by M. larson and F. Bengzon,  and on this book to find a posteriori error estimate for this equation by linear fem, they used jump for $\frac {\partial u}{\partial n}$
$$\eta_k=h_k\| f+\Delta u_h\|_{L^2(k)}+\frac{1}{2}h_k^{\frac{1}{2}}\|[n.\nabla u_h]\|_{L^2({\partial k}-\partial{\Omega})}$$

Title: Jump conditions on interior boundaries


问:
String search algorithm on small alphabet
I am looking for an algorithm to search a substring in a string.  I know there are a number of wellknown Algorithms like Knuth-Morris-Pratt, for instance, and I suppose most preimplemented  functions use one of those. However, I remember a lecture, long ago, where the lecturer said that the usefullness of such an algorithm depends partly on the size of the alphabet. KMP and BMH obviously work fine on an ordinary 26-letter alphabet. But what do I do if I have a much smaller alphabet (say, DNA: 4 letters or simply a binary)?

Are there an particular algorithms that work well on very small alphabets?
I googled and supposed I should easily find something as searching through DNA is rather common, but to no avail. Any help?

Please do not downvote, this is my first question in this subcommunity, I do not know yet what is considered too basic here and what is ok.


问:
Periodic boundary conditions for solving Navier Stokes Equations on a Staggered Grid
I want to solve two dimensional Navier Stokes equations on a staggered grid for the case of Taylor-Green Vortex. My initial conditions are standard sine and cosine functions. As I am aware, I should use periodic boundary conditions on all the four boundaries. I am facing some trouble regarding implementation of this BC as I could not find enough literature about it.

Suppose, if my index ranges from $0$ to $N$. Then, is this the correct way to use periodic BC?

$$u_\text{vel}(N) = u_\text{vel}(1)\\
u_\text{vel}(0) = u_\text{vel}(N-1)$$

and similarly for $v_\text{vel}$ and pressure variables?


问:
Conjugate gradient: the 1-norm of the residual
I am trying to solve $Ax=b$ using the conjugate gradient method. However, it is important to me to obtain a bound not only on the usual residual $||b-Ax_k||_2$ but also on the quantity $||b-Ax_k||_1$.

I have two questions:

Is it true that $||b-Ax_k||_1 \leq C ||b-Ax_0||_1$ for all iterations $k$ where $C$ is some constant?

In other words, is it true that the conjugate gradient amplifies $1$-norm of the residual by at most a constant factor?

If we can choose $C=\sqrt{n}$ (where $A \in \mathbb{R}^{n \times n}$) then the above fact is trivial since its true with the 1-norm replaced by the 2-norm and using the bound $||\cdot||_1 \leq \sqrt{n} ||\cdot||_2$. However, I'd like to $C$ to be a constant not depending on dimension $n$.

If the answer to question one is negative, is there some modification of the conjugate gradient which achieves such a bound?


问:
Implementing std::nextafter: Should denormals-are-zero mode affect it? If so, how?
This might be the wrong stackexchange site for this question.  math.SE, cs.SE, programmers.SE, and of course stackoverflow are all possibilities.  I'm hoping to reach an audience that might actually use this function, to get some feedback on how they'd like it to work.

C's double nextafter(double x, double y) function returns the next representable value from x, in the direction of y.  (commonly + or - infinity, but can be anything).

I'm tidying up GNU libc/libm's implementation to compile to nicer code.  (It currently extracts each double to two 32bit integers.  It's pretty clunky and uses lots of branches, and doesn't take advantage of int64_t.).  If I'm going to change it at all, I might as well make as much improvement as possible, but using an FP compare-equal would change the results in DAZ mode.

How should it behave when the denormals-are-zero and/or flush-to-zero are active?  On x86, there are two separate flags which sacrifice gradual underflow for performance.  By default, neither is set.  One or both can be set (they're often used together, but don't have to be):

DAZ is an input filter: When an FP math / compare instruction reads its inputs, denormals are considered +/-0.0.  So a compare between two denormals finds they're equal.  Arithmetic can easily produce denormal results, though.

FTZ is an output filter: Denormal inputs work normally, with no effect on compare instructions.  However, ordinary math instructions can't generate denormals.  Denormal results are flushed to +/-0.0.  DBL_MIN/2 + DBL_MIN/1.5 should give the same result as without FTZ (e.g. if those input denormals are loaded as constants).

In both cases (esp DAZ), one might argue that the next representable value below DBL_MIN (the smallest normal number) was +0.0, and that the next number below +0.0 or -0.0 was -DBL_MIN.  Implementing this would require checking the FTZ/DAZ flags.

One might argue instead that in DAZ mode, two different denormals should be considered equal, because they do compare equal.  Thus, it's correct to return y when x equals y, as required by POSIX, and by the C11 standard (7.12.11.3 pg 256 and Annex F.10.8.3 pg 529).  (The motivation for this specification is for nextafter(0.0, -0.0) to return -0.0, and vice versa.)

The 2nd behaviour allows an extremely fast & compact implementation (no taken branches and 3 not-taken branches polluting the branch-predictor, and ~9 clock cycles latency on Intel Haswell (for comparison, an FP multiply takes 5 cycles) for the common case (x!=y, and x!=0, and the result doesn't overflow(infinity) or underflow (denormal)).  Less than half the total code size compared to glibc's current version.

The C11 standard says:

"Even though underflow or overflow can occur, the returned value is independent of the
  current rounding direction mode."

But DAZ isn't a rounding direction mode.

The boost::math::ulp function's documentation observes that "our experience is that std::nextafter often breaks depending which optimizations and hardware flags are in effect".

Consistency with other implementations is desirable for glibc's nextafter function, but consistency with glibc's existing behaviour may take precedence.  The current implementation is not influenced at all by any FP settings, because it uses only integer comparisons on the IEEE float bit patterns.  Consistency with nexttoward, and with nextafterf(float,float) / nextafterl(long double, long double) is also important.  I haven't looked into the long double versions yet.  Integer compares are clunkier there, because int64_t can't hold 80 bits.

Using integer-only compares is only slightly slower / bulkier to implement (using a 64bit integer compare, and the sign of x and y, instead of just an FP compare and the sign of x).  IEEE floats are cleverly designed so their bits can be compared as sign-magnitude integers (and thus as 2's complement integers with a check for the both-negative case).

How does nextafter behave on other systems (OS X, MSVC++, anything else)?  If people are interested in testing it, I might write a little test program.

How do you wish nextafter behaved?  (nextafter is pretty much fully specified except for DAZ/FTZ, so wishes other than that might be better fulfilled by different functions.)

Did I miss anything?

I know the current implementation is "good enough" for such a rarely-used function.  It's ugly, but it's used so rarely that it doesn't have much impact.  (Actually, is it rarely used?  Has anyone ever seen code that used it a lot?)

Part of the motivation here is that I enjoy optimizing things in assembly language, and integer ops on FP data is interesting, especially on x86-64 where integer instructions can be used on the vector registers that are also used for floats and doubles (but not 80bit long double for nextafterl, or for the direction arg of nexttowards).


问:
Evaluation of feynman diagrams: FORM vs. Sympy vs. other
I just started working on my master thesis in the theoretical particle physics and my supervisor asked me to start working with FORM, a script interpreter with basic knowledge of lorentz and dirac structures implemented, for the evaluation of some feynman diagrams. As my first glimpse on how FORM looks and works it doesn't seem to be very user friendly and intuitive. So I started to look around a bit and found that sympy, which I know the basics of, also has a module for lorentz indices and dirac matrices.

So before I work more into one of these I thought I would seek some solid advice from someone who maybe has worked with both or was once confronted with the same decision. Has there ever been a benchmark of these and if FORM is the faster option is it worth the inconveniences? Is there maybe a third option considered superior over the two mentioned?

I thought that would be on-topic here, but please correct me if I'm wrong.


问:
Optimal distribution of zeros and ones over matrix
I have the following problem:

Given a matrix with n rows and m columns. Some elements of the matrix are unavailable. For each column, you have a set containing a number of zeros and ones which must be distributed over the available elements in that column. The total number of zeros and ones in a set is equal to the number of available elements in the corresponding column. The question is: How should the values of the sets be distributed over the elements in the corresponding colum, such that the number of rows for which the total number of ones is less than 20% of the available elements in that row, is minimized.

Does someone knows a similar problem which has a polynomial time algorithm solving the problem to optimality, or knows a polynomial time algorithm solving the problem to optimality (if it exists)?


问:
Stability in Discretization of 1D Stationary Boltzmann equation
I want to discretize and numerically solve the following PDE:
\begin{equation}
v(k)\dfrac{\partial f}{\partial x} + E(x)\dfrac{\partial f}{\partial k} = S\{f\}
\end{equation}
using finite volume (box integration) method.
Since the direct discretization results in an oscillatory solution, I have used the upwind discretization scheme. This eliminates the oscillations and the solution satisfies the maximum principle. However, if I have my $x$-grid as different $x_i$ points with $i=1,2,...,N$, the upwind discretization schemes loses the connection between $f(x_1)$ (the boundary value) and $f(x_2)$ because I am always calculating the right-sided difference on the second grid node.

When I try to use the central difference on $x_2$, oscillations and instabilities return to the solution, and when I use the upwind scheme, the boundary points are disconnected from their adjacent nodes.

Any suggestions on how to solve this issue?


问:
Admissible box constraint for a quadratically constrained linear program
I am looking at a real-world resource allocation problem that is cast as a quadratically-constrained linear program of the form

$$
\max\langle f,x\rangle
$$
subject to
$$
\begin{aligned}
m \leq\,\, & x \leq M, \\
A\cdot x & \leq b, \\
A_{eq}\cdot x & = b_{eq}, \qquad\text{and} \\
\langle x,Q\cdot x\rangle & \leq \alpha.
\end{aligned}
$$

It is a small problem, that is, $x\in\mathbb{R}^{55}$, and it can be solved in much less than a second using Matlab's $\verb|fmincon|$ to obtain an optimal solution $x^*$. The question is the following: is there an "intelligent" (read: systematic/efficient; again, hacks are welcome) way to determine an admissible box constraint of the form
$$
\tilde{m}\leq x*\leq \tilde{M}\ ? \qquad(1)
$$
That is, if $B$ is the region described by the box constraint (1) and $F$ is the feasible region of the original problem, then it should be the case that
$$
B\subseteq F.
$$
It is fine if $B$ has empty interior, as long as the relative interior with respect to the affine subspace that contains $B$ is not empty.

To provide some context, the components of the optimal solution $x^*$ are communicated to production departments that operate independently and $B$ represents allowable deviations from their targets. Let me point out that, as is implicit in the description of the problem, the original constraints are hard constraints and it is not meaningful to violate them; i.e., to have a box $B$ around $x^*$ that extends beyond the feasible region.

Question: Any pointers to the literature, high-level hints to the right direction, or quick and dirty approaches to the following problem will be highly appreciated.


问:
The numerical solution of a (very ugly) set of integro-diferential equations
I've been having fun playing around with a simple model for reaction in a co-flow fluidized bed (I'm afraid I'm a chemical engineer, not a computer scientist). Unfortunately, simple physical assumptions often lead to terrifying equations. I've reduced the problem to solving the following set of simultaneous integro-differential equations for the 2 unknown functions $g(s)$ (on $s\in[s_0,s_{max}]$) and $C(z)$ (on $z\in [0,L]$), where $f(C,s)$ is a known function and $\alpha$ and $\tau$ are constants:

$$-\frac{dC}{dz} = \alpha \int_{s_0}^{s_{max}} f(C(z),s) g(s) ds$$
$$g(s) = \frac{\exp\left(-\frac{L}{\tau}\int_{s_0}^{s} \frac{ds}{\int_{0}^{L}f(C(z),s)dz}\right)}{\frac{\tau}{L} \int_{0}^{L}f(C(z),s)dz}$$
and the one BC:
$$C|_{z = 0} = C_0$$

So far I've tried the brute-force approach: finite differences. I divided the intervals  $[0,L]$ and $[s_0,s_{max}]$ into $n$ points, and defined the $2n$ unknown variables $C_1, C_2,\ldots,C_n$ and $g_1, g_2,\ldots,g_n$, which are just the values of my unknown functions at these points. I then approximately the derivatives and integrals above using FD expressions, and attempted to solve the system of $2n$ non-linear algebraic equations using $\texttt{fsolve}$ in MATLAB.

Perhaps predictably, MATLAB doesn't like this, and fsolve struggles to give a solution for medium-sized $n$ ($n=100$ returns 'no solution found', and $n = 250$ is still running after half an hour.) My questions would be:

Are there any ways of solving large systems of simultaneous non-linear equations, that are more efficient than $\texttt{fsolve}$. Preferably I'd be looking for an algorithm already implemented in MATLAB.

Even more optimistically, are there any other possible approaches to solving these types of highly non-linear integro-differential equations?

Thanks


问:
Generating Random Numbers in Fortran for Metropolis method
I am having a really hard time getting any kind of reliable / consistent result from my Metropolis code. I have torn it apart and am now examining just the randomness in my random number generator.

I decided I would just run many iterations over a spread from -1 to 1 and it should give me an average of 0. I was completely surprised by the lack of precision. My code does give results about zero, but, from -200 to 200 with almost an equal distribution throughout the range.

I'm trying to find a way to bring this closer together. I have tried more runs (knowing 10,000 isn't really that many) but even up to a billion runs, the averages can range just as much if not more.

Any suggestions on how to fix this randomness?

I have been using:

```<code>    program test_rand
    real*8 x, y, z, EAx, EAy, EAz, Tx, Ty, Tz
    integer j, i

    CALL init_random_seed

    Do i = 1, 10000
        EAx = 0.D0
        EAy = 0.D0
        EAz = 0.D0
        Tx = 0.D0
        Ty = 0.D0
        Tz = 0.D0
        Do j=1,10000
            CALL RANDOM_NUMBER(x)
            CALL RANDOM_NUMBER(y)
            CALL RANDOM_NUMBER(z)
            EAx = 2.0D0 * x - 1.0D0
            EAy = 2.0D0 * y - 1.0D0
            EAz = 2.0D0 * z - 1.0D0
            Tx = EAx + Tx
            Ty = EAy + Ty
            Tz = EAz + Tz
        end do
        write(50,*) i, Tx, Ty, Tz
    end do
    end program test_rand

    !   initialize a random seed from the system clock at every run (fortran 95 code)
    subroutine init_random_seed()
    INTEGER :: i, n, clock
    INTEGER, DIMENSION(:), ALLOCATABLE :: seed
    CALL RANDOM_SEED(size = n)
    ALLOCATE(seed(n))
    CALL SYSTEM_CLOCK(COUNT=clock)
    seed = clock + 37 * (/ (i - 1, i = 1, n) /)
    CALL RANDOM_SEED(PUT = seed)
    DEALLOCATE(seed)
    end
</code>```


问:
Are there any benefits of computable analysis to numerical algorithms
Computers can work only with computable numbers, while most of the algorithms are based on analysis of real numbers (real analysis).

When I heard of the existence of computable analysis I immediately had these questions:

Can it help numerical analysis in some way?

Can we for example, help prove convergence of algorithms that are not convergent in the reals?

Can we find different algorithms that work for computables faster than they do for reals?


问:
Calculus of Variations with unknown cost function but some data
I have a problem that I've framed out in a particular way, but I don't know if I'm re-inventing the wheel here. Is there an existing literature base in this problem? Does it have a corresponding term of art? Where would you recommend that I start reading about this type of problem?

The basic problem that I have is that I need to make predictions about where ants will be in a large terrain. I have all sorts of information about the terrain -- its elevation, bodies of water, the flora cover, the distances to cities and other human settlements, the weather -- all of these things plausibly influence where and how ants move across this terrain.

In addition to the terrain data, I have data about where ants have historically been found with lat/longs to a high degree of precision.

We also know that ants will generally be moving from one end of the terrain to the other. But there is not a hard requirement that they always move closer to the opposite end at every step. We just know that they're starting at one end and want to go somewhere on the other end. The "box" of terrain is just a very wide rectangle, and I'm comfortable with assuming that ants want to go from the longest side to the opposite longest side.

Plausibly some ant could go from one corner to the farthest corner, but it's more likely that they'll go to a nearer point because the path is shorter, so probably the path is less costly. This is where the calculus of variations comes in: I want to make inferences about the plausible paths that ants will take. Likewise, an ant could make a several U-turns to avoid some obstacle (like a suddenly steep elevation, or crossing a body of water), and doing so would incur a cost cost, but probably won't make several such U-turns, if judicious planning could avoid it ahead of time by taking an alternative route. This is similar to how water will usually flow down hill, but special circumstances can make a geyser.

Conceptually, the cost function is a function of the terrain data, but the cost itself is unknown.

The straight-up machine learning approach would just be to use the terrain data to make inferences about where ants will be directly as some sort of regression problem -- this can be done any number of arbitrary ways, but I wonder if doing so without considering paths will lead to implausible inferences, such as making predictions that ants might be at any point along a line, except some vast swath in the middle where the model predicts probability 0, or that there is a "dotted line" of alternating high- and low-probability regions. I wouldn't really believe that model, and would tend to think that it's overfitting the data. It seems like making inferences about the cost function, however indirectly, would impose a certain amount of regularization on the solution driven purely by the "pure" machine learning regression problem, so incorporating that information would improve our inferences.

So where would I start reading about this style of problem?


问:
Optimized open source BLAS / LAPACK package
I was wondering what is a more optimized open source BLAS/LAPACK package with respect to modern multi-core processors (Haswell and beyond).  Is there any distribution that can attain performance close to that of Intel MKL for instance?

In addition to Intel, which libraries may perform well on AMD architectures for Zen gen1 and Zen gen2 ?

I was hoping that the answer can be based on actual performance analysis work.

Thanks
Michael

ContentLicense: CC BY-SA 4.0


问:
Numerical Euler Rotation Equation
The problem I have may be really simple, but still getting a hard time solving it. So I have the Euler rotation equations:

$$I_{1}\dot{\omega}_{1}+\left(I_{3}-I_{2}\right)\omega_{2}\omega_{3}=\lambda_{1}$$

$$I_{2}\dot{\omega}_{2}+\left(I_{1}-I_{3}\right)\omega_{3}\omega_{1}=\lambda_{2}$$

$$I_{3}\dot{\omega}_{3}+\left(I_{2}-I_{1}\right)\omega_{1}\omega_{2}=\lambda_{3}$$

where the $I_{i}$ are the moments of inertia about the principal axes of rotation and $\omega_{i}$ the time dependent angular velocities about each axis. In general $I_{i}\neq I_{j}, i\neq j$, or don't satisfy any situations in which they can be reduced to easier relations. $\lambda_{i}$ is a normal distributed random number.

I know that these equations are non-linear and, in general, have no analytic solution. The question is:

$\qquad$ Is there a good numerical integrator for these equations?

I'm using C++ and have looked into the LAPACK package, but I'm kind of confused on how to use it.

I know the reference from Skowron and Gould (arXiv:1203.1034 [astro-ph.EP]); however I really don't know how to implement this algorithm.

If somebody could help me in finding an open source integrator or a reference where they talk about the implementation of this code, it would be great.

Title: Numerical Euler Rotation Equation


问:
Representing a 3D system in 2D (Electromagnetic modelling)
Ok so I'm a complete beginner in computational modelling (I use analytical methods of physics typically) but I would like to model an anisotropic, aperiodic (but not random) finite array of metallic nanorods, all of which have their long axis orientated in the same direction as eachother, but this direction can be completely arbitrary.

What I am interested in is the charge oscillations in each nanorod (localised surface plasmons) which can electromagnetically couple to neighbouring nanorods and generate a collective charge density wave (plasmon/polariton) which can in principle propagate through the system. For example I could impose a periodic driving of the charge density in one nanorod at a particular frequency and see how it decomposes into the different allowed momentum states of the system and propagates. I want to have a simulation where I can see the propagation of these waves.

The problem is that because of the aperiodicity I can't use periodic boundary conditions. I believe this means that a full 3D electromagnetic simulation is too computationally heavy (I have a good spec iMac but nothing better). So I assume I have to do a 2D model, but I have no idea if this makes any sense for my system and how to implement. These nanorods in principle can have any orientation in 3-dimensions and the EM field is intrinsically 3D so how can I map this problem to 2D. How reliable is this approach for modelling the reality?

I'm hoping for some suggestions on what software to use and how to get started and perhaps some insight into this 3D to 2D mapping problem. I have access to COMSOl but happy to try other (free) software.


问:
Stability of two PDEs
I have the following two PDEs which I want to check for stability:

$$u_t= u_{xx} , \ u(x,0)=1 , \ u_x(1,t)=-hu^4(1,t) , \ u_x(0,t) = 0 $$

$$u_t = u_{xx}-\sin(x+t)+\cos(x+t) , \ u(x,0)=\cos(x) ,\ u_x(1,t)=-hu^4(1,t)-\sin(1+t)+h\cos^4(1+t) , \ u_x(0,t)=-\sin(t)$$

For the first I got that $d/dt (u,u) = d/dt \int_0^1 u^2 dx = 2[-hu^5(1,t)-\int_0^1 u_x^2dx$, but how to estimate $u^5(1,t)$?

Just before I forget, h is some parameter greater than zero.

Thanks!

Title: Stability of two PDEs


问:
Calculating lagrange polynomial for 100 points?
I need to calculate the lagrange polynomial  which approximates $e^x$ at $101$ points, the points $\frac{k}{101^2}$ for $k\in\{0,1,2\dots 100\}$.

I tried the following code:

```<code>import java.math.*;
import org.apache.commons.math3.analysis.polynomials.PolynomialFunctionLagrangeForm;

public class interpolation {

    public static void main(String[] args) {
        // TODO Auto-generated method stub
        int N = 100;

        double [] X = new double [N+1];
        double [] Y = new double [N+1];

        // Here we give the computer the points:

        for(int k=0; k<=N ;k++){
            double a = k/(N*N*1.0);
            double b = Math.exp(a);
            X[k]=a;
            Y[k]=b;
        }

        // here we build the polynomial

         PolynomialFunctionLagrangeForm pol = new PolynomialFunctionLagrangeForm( X , Y );

         double [] P = pol.getCoefficients();

         // here we print the polynomial

         for(int i=P.length -1 ;i >= 0;i--){
             System.out.print( P[i] + "    ");
         }

    }

}
</code>```

It gives good results for $N=5,10$ but bad ones for $N=50,100$ (it doesn't resemble the Taylor expansion at all). I would like help in calculating this polynomial more precisely. Thank you kindly.


问:
Is there a difference between the Galerkin MWR and other techniques called Finite Element Method
Recently, I took a course on numerical methods where we learned about Galerkin's Method of Weighted Residuals. We were told that it forms the basis of the Finite Element Method .

However, in most introductory texts I've seen on FEM, they start with a simple example where they manually write out the equations to minimize the potential energy at each node. For instance, see this example of a system of springs. There are no test functions, so it doesn't seem to be the MWR.

So, my question is:

Is there a name for this other method?

What are the practical and numerical differences between it and the Galerkin Method of Weighted Residuals? (e.g., is it easier to set up for certain problems? Does it have better numerical precision or stability?)


问:
Alternative to two "for" loops in finding best neighborhoods for TSP?
I am trying to solve Travelling Salesman Problems using tabu search. I have been able to successfully find "near enough" optimal solutions (as well as one optimal, yay!).

For the moment I am using sub-tour reversal as the local search precedure. The basic gist of it is similar to 2-opt, but rather than swapping two cities, it reverses the subtour between and including those two cities. I.e if we have a tour [1, 2, 3, 4, 5, 6, 1] and we swap city 2 and 5 we get [1, 5, 4, 3, 2, 6, 1].

However I am using two "for" loops to find the best neighborhood from the initial tour. But as you know, going through every single swap is a very time consuming process, which limits me to TSP problems of less than 500 cities. Anything above 1000 cities will take my program hours.

So my question is, can I apply some other method than two or loops to find the best neighborhood in more efficient (fast)?

For now my "for" loops look like this:

```<code>int col = 1;
for (int i = 1; i < initialSolution.length-1; i++) {
    for (int j = col; j < initialSolution.length-1; j++) {
        if (i == j) {
            continue;
        }
        if ((i == 1) && (j == initialSolution.length-2)) {
            continue;
        }
        if ((j == 1) && (i == initialSolution.length-2)) {
            continue;
        }
        // LOCAL SEARCH PRECEDURE
    }
}
</code>```

such that I do not repeat previously reversed tours, and such that I do not reverse the whole tour.


问:
Gradient convergence on a checkerboard domain in finite element
I've got an elliptic problem in a domain with two materials in a checkerboard arrangement. Both materials are highly dissimilar, being one much more compliant than the other.

If I use classic linear lagrange finite element, because of the corners of the checkerboard, there will be singularities that will make the solution gradient not converge with the mesh refinement. How can I prove this fact? I am looking for a proof based on the fact that the material properties uses the same mesh discretization, but it's piecewise constant. Which references would be good for this proof? Thanks.


问:
Nonlinear conjugate gradient restart threshold 1/10
Nocedal and Wright
on Conjugate Gradient Methods, p. 123, describe a

restart strategy ... whenever two consecutive gradients are far from orthogonal

$\qquad {{| \nabla f_k^T \ \nabla f_{k-1} |}
\over { \|\nabla f_k\|^2 }} \ge \nu $,
with $\nu$ typically 1/10.

Can anyone comment on CG with such restarts, or point to test cases on the web ?

Or is the "popular choice max $( 0, \beta^{PR} )$"
(Nonlinear_conjugate_gradient_method)
good enough, satisficing ?

(A good answer to bfgs-vs-conjugate-gradient-method says,

Anecdotal evidence points to restarting being a tricky issue,

  as it is sometimes unnecessary and sometimes very necessary.

Well, that's generally true of a lot of things (taxes come to mind).

Test cases with plots of $\beta_k$ or $\theta_k$ might be interesting.)

A possibly silly test case that led to the question is CG on an ill-conditioned quadratic in 2d:

```<code>import numpy as np
from scipy.optimize import fmin_cg

n = 2
cond = 100
eigenvalues = np.linspace( 1./cond, 1, n )
xmin = 1000 * np.ones( n )

def fprime( x ):
    return eigenvalues * (x - xmin)
def f( x ):
    return (x - xmin) .dot( eigenvalues * (x - xmin)) / 2

x0 = np.zeros( n)

ret = fmin_cg( func, x0, fprime )
</code>```

Added:

for n in [1,2,3,4,5] this takes 80 81 6 40 9 iterations,
721 722 28 94 27 function evaluations.
(Is CG generally very sensitive to the linesearch ?)

The Mathematica conjugate gradient minimizer
has a RestartThreshold with default 1/10.
But, sorry, I don't speak Mathematica.
Any native speaker care to try this ?


问:
PETSc - difference between PetscNew and PetscMalloc
I wonder what is the difference between both functions.

In C++, the difference between new and malloc is easier: new allocates memory by calling the constructor.

What about PETSc ? In the example, it seems to be pretty close:

in dm/examples/tutorials/ex6.c:

```<code>struct _p_FA {
  MPI_Comm   comm[3];
  PetscInt   xl[3],yl[3],ml[3],nl[3];    /* corners and sizes of local vector in DMDA */
  PetscInt   xg[3],yg[3],mg[3],ng[3];    /* corners and sizes of global vector in DMDA */
  PetscInt   offl[3],offg[3];            /* offset in local and global vector of region 1, 2 and 3 portions */
  Vec        g,l;
  VecScatter vscat;
  PetscInt   p1,p2,r1,r2,r1g,r2g,sw;
};
typedef struct _p_FA *FA;
FA         fa;
PetscNew(&fa);
</code>```

In  ksp/ksp/examples/tutorials/ex42.c:

```<code>PetscMalloc(sizeof(struct _p_CellProperties),&cells);
</code>```

The rule is just: using PetscMalloc for a structure, and PetscNew for a pointer of structure?


问:
Preconditioned Steepest Descent
For my program assignment I need to write a preconditioned steepest descent algorithm. I have psedo-code from my professor which is here:

From reading in my text on this method they say that $P$ is the preconditioned matrix but they do not define it in any way shape or form. I am not sure at all, but by my professors psedo-code I would guess that $M$ is $P$? I honestly have no idea if that is the case. I have tried to search online for this but I seem to be getting information that is not very consistent with the way our book/professor presents the material.

Update: $M$ is just a general precondition matrix; i.e. Jabobi, or any other precondition matrix. The problem I seem to be having now is that I do not really understand the pseudo-code. I have all the ingredients necessary to start but I am not sure how to begin.

Note: the book I am using is,  A. Quarteroni, R. Sacco, and F. Saleri, Numerical Mathematics,

Title: Preconditioned Steepest Descent


问:
CFD implementation in software
I've been learning CFD theory for 4 months (which includes FEA, FDM AND FVM) and now I want to start running simulations. As a novice to CFD software, I would appreciate some advice on where to start, e.g., OpenFoam, Matlab, Fortran? (I am comfortable with programming and learning coding languages). </p>

Title: CFD implementation in software


问:
Exact line Search in Steepest descent
I wanted to clarify the idea of the exact line search in steepest descent method.

An exact line search involves starting with a relatively large step size ($\alpha$) for movement along the search direction $(d)$ and iteratively shrinking the step size until a decrease of the objective function is observed. Is this correct?
I.e. the only condition suppose to be satisfied is: $f(x_{k+1})<f(x_k)$, where $x_{k+1} = x_k + \alpha_kd_k$

In exact line search the step size $\alpha$ is not necessary to be fixed. Is this true?

If this is true then what is the difference of the exact line search from the backtracking line search? I know that backtracking line search is based on Armijo-Gooldstein condition or Wolfe's condition.

Title: Exact line Search in Steepest descent


问:
Check the well-posedness of the problem with level set equation
Given a domain $\Omega$ with 2 subdomains $\Omega_1,\Omega_2$ as figure below

Level set equation $$ \partial_t\varphi + U\cdot\nabla\varphi =0 $$

$U$ is velocity only presents inside $\Omega_1$. We assume $U$ is irrotational, i.e. $\nabla \times U=0$, which implies that the velocity ﬁeld can be derived from a potential function $\Phi$ defined in $\Omega$ such that $U=\nabla\Phi$. The potential equation is given by

$$
\begin{cases}
-\nabla\cdot(\nabla\Phi)&=-\beta_S \quad \text{in }\Omega_1\\
-\nabla\cdot(\nabla\Phi)&=0 \quad \text{in }\Omega_2
\end{cases}
$$ with $\beta_S$ given.

I want to solve this level set equation, especially with unfitted mesh but I don't know the problem is well-posed or not.


问:
Best way to add a positivity constraint to Newton's Method
So given an objective function $f({\bf x})$, I would like to include a positivity constraint when I perform the fixed point iteration:
$${\bf x}^{(t+1)}={\bf x}^{(t)} - \text{H}_f^{-1}\nabla f({\bf x}^{(t)}),$$

for ${\bf x}^{(0)}>{\bf 0}$.

The best way I've found to do it is to add a log barrier to the objective function so it becomes
$$g({\bf x}) = f({\bf x})+\frac{1}{n}\sum_{i=1}^n\log{\bf x}_i,$$
since this will cause $\nabla g$ to shoot to infinity as one approaches $\bf 0$ from the positive orthant.  And then to modify the iteration step size to prevent the step size from being large enough to cross all the way over the log barrier to the other side:
$${\bf x}^{(t+1)}={\bf x}^{(t)} - \frac{1}{2}\min\Big[2, \frac{\left\|{\bf x}^{(t)}\right\|}{\left\|\text{H}_g^{-1}\nabla g({\bf x}^{(t)})\right\|}\Big]\text{H}_g^{-1}\nabla g({\bf x}^{(t)}).$$

Is this the best way to enforce a positivity constraint on Newton's Method, or is there a better way?


问:
Scaling a vector-valued non-linear function for numerical optimization/minimization
I am trying to minimize a non-linear vector-valued function in MATLAB.  As a test case for my code, I try to minimize a function whose solution I know apriori.

The problem is that one of the solutions is widely distributed in scale.  They are [175,164,854,3.7e5,6000];

As you can see, the 3.7e5 has a wide difference in scale with the rest of the solution vector. As a result, despite changing tolerances, algorithms etc., the best that MATLAB can find for my x(4) solution is 5.2e5

Now, I know that using a diagonal (or close to diagonal) scaling can help in finding the minimizing solution. However, from my understanding, the Hessian matrix is needed to obtain the scaling factor, i.e. Ideally I would like to scale the solutions so that the Hessian of the objective function becomes the identity matrix at/near the solution.

But the problem is, the Hessian of a vector-valued function is a Tensor and not a matrix.

Is there a solution, or any other method for handling this problem ?


问:
Solving a nonlinear poisson equation via variational minimization
I am kind of new in finite elements and I am solving simple "Poisson nonlinear" problem.

$- \nabla ((1 + u^2) \nabla u) = f$

$u = 0 \ \text{on} \  \Omega $

I am using Newton solver, where I have defined step as following

$J = s * ( -g )$,

where  $ J = \int (2 \delta u \nabla u_h) \nabla v + (( 1 + u_h^2) \nabla \delta u ) \nabla v  dx $

and

g = $\int (1 + u^2 ) \nabla u \nabla v - f v \ dx$.

I would like to check my solution and monitor energy, however I am bit stuck here. I know that in linear case it would be QP function, but how is my strain energy function defined here?


问:
SPECT reconstrction using MLEM
In  Single-Photon Emission Computerized Tomography (SPECT) parallel beam reconstruction using Maximum-Likelihood Expectation–Maximization(MLEM), is it sufficient to scan the object around 180 degree?  or it is required to perform full 360 scan?</p>

Title: SPECT reconstrction using MLEM


问:
What is the best option in terms of library or software to solve this system of hyperbolic PDEs?
I want to solve a system of coupled nonlinear 1-D PDE $(\partial_{tt} + \alpha\partial_t)u_i(x,t)=\partial_{xx}(\sum_{j=1}^{j<i}ju_j(x,t)+i\sum_{j=i}^{n}u_j(x,t))-\sin(u_i(x,t))+f$, using method of lines and for this I defined a spatial grid of about $N\approx1000$ points. Since my PDE is second order both in time and space, the number of ODEs that I have to solve is about $2n*N$ which needs to be run for enough time to get the solutions. For this large system of coupled ODEs I used an ODE solver (tested both implicit and explicit) and although the fastest is ODE23, yet it takes a long time when my $\delta x=0.01$. Since I have access to a multicore system, I was wondering if I could solve this system in a scalable method. I do not think it is applicable using MATLAB or MATHEMATICA. What other options do I have?</p>

```<code>% spatial grid
x=xl:dx:xu;
%intial condition
u0(1:n) = somdefunction (x)
u0(n+1:2n)=0;
%ode solver
[t,u]=ode113(@fun,tspan,u0,options)
% pde
function ut = fun(t,u)
uxx = 4*del2(u(1:n),dx);
%boundary conditions
% pde
ut1(1:n) = u(n+1:2n);
ut2(1:n) = uxx - sin(u);
ut = (ut1;ut2);
end
</code>```


问:
Galerkin FEM error when using even number of elements

Intro:

I am developing a Galerkin FEM code in matlab - starting small with a simple 1D ODE. The equation I'm trying to solve is: $$a u_x = cos(x), x \in [0, 2\pi]$$ Which has a known exact solution $$ u = sin(x)/a $$

When I solve the equation using my code on an odd number of linear elements (nodes = N+1), the solution is as expected:

Problem

I'm running into a problem when I set the number of elements to an even number, in which case the solution vector is just filled with NaNs and I receive the following warning/error:

Warning: Matrix is singular to working precision

Most of the questions I've looked up on mathworks answers have to do with incorrectly constructed matrices - that's not the case here. I'm not sure if I'm running into a numerical issue, where a zero is too small, or if this is an implementation issue.

I construct both the elemental stiffness matrices and the RHS vector in the same way, and I can't tell a difference between the two systems from a quick comparison of their respective Ax=b structures.

Methodology

I'm generating the elemental stiffness matrices individually, where N is a basis function. The weak form of the ODE above for an elemental stiffness matrix is:
$$ I_e = \int_\Omega \underline{N} a u_x dx = \int_\Omega \underline{N} cos(x) dx $$

Using integration by parts (i.e. magic), this equation turns into this:
$$ (au) \underline{N} |_{\partial \Omega L}^{\partial \Omega R} - \int_\Omega (au) \underline{N}_x dx $$

Then I substitute a continuous value of u for the weighted sum of node values ($\hat{u}$), as such:
$$ u = N_1 \hat{u}_1 + N_2 \hat{u}_2 = \underline{N}^\intercal \underline{\hat{u}} $$

Which is finally equal to the following:

$$ I_e = a \underline{N}^\intercal \underline{\hat{u}} \underline{N} |_{\partial \Omega L}^{\partial \Omega R} - a \int_\Omega \underline{N}^\intercal \underline{\hat{u}} \underline{N}_x dx  $$

Since the first term in $I_e$ cancels out except for the boundaries, in which case $u$ is already 0, I'm assuming that term can be safely ignored. The resulting elemental stiffness takes the form:

$$ I_e = \begin{bmatrix} -a \int_\Omega N_1 N_{1x} dx & -a \int_\Omega N_1 N_{2x} dx \\ -a \int_\Omega N_2 N_{1x} dx & -a \int_\Omega N_2 N_{2x} dx \end{bmatrix} $$

Code

The entire code is in a repository on my github here, but the most important part is here below. All of the functions called from this script have been checked and verified to work properly, namely basis_1D(), lgquad(), and lgwt(). Feel free to check the repo if you're curious

The portion of the code that constructs A and b is:

```<code>for kk = 1:numElem
    % Initialize elemental stiffness matrix
    Ie  = zeros(2);

    % Anonymous function that maps natural coordinate to global coordinate
    X   = @(s)basis_1D(1,1,s,0)*coordPairs(kk,1) + ...
                basis_1D(1,2,s,0)*coordPairs(kk,2); % X = x(s)

    % Anonymous function that calculates the Jacobian as a function of the natural coordinate
    J   = @(s)basis_1D(1,1,s,1)*coordPairs(kk,1) + ...
                basis_1D(1,2,s,1)*coordPairs(kk,2); % J = dx/ds

    for ii = 1:2
        for jj = 1:2
            % Anonymous function to what I denote above as N_x
            fun1    = @(s) basis_1D(1,ii,s,1)./J(s);

            % Anonymous function to what I denote above as N
            fun2    = @(s) basis_1D(1,jj,s,0);

            % The two functions combined together, as well as the jacobian for numerical quadrature
            fun     = @(s) fun1(s)*fun2(s).*J(s);

            Ie(ii,jj) = -a*lgquad(fun,[],-1,1); % Calling Guass-Legendre Quad routine

            II = (kk-1)+ii;
            JJ = (kk-1)+jj;
            A(II,JJ) = A(II,JJ) + Ie(ii,jj); % Assembling global matrix
        end

        fun1    = @(s) cos(X(s));
        fun2    = @(s) basis_1D(1,ii,s,0);
        fun     = @(s) fun1(s)*fun2(s).*J(s); % Combined function for RHS

        % Calling Gauss-Legendre Quad routine and assembling RHS
        b((kk-1)+ii) = b((kk-1)+ii) + lgquad(fun,[],-1,1);
    end
end
</code>```

Can someone give me an idea of what my issue might be?


问:
Convolutional Neural Nets for Large Scale Image Classification on Windows?
I need to build some CNN's using very large numbers of photos, to classify which photos contain certain features.

I've run into a number of road blocks.

As much as I love both Linux and virtualization, I've found that Linux-only solutions like Berkeley's Caffe and NYU's pre-trained OverFeat don't work for large scale data aren't viable because virtualization has a lot of problems in an HPC and big data scenario, especially when the tools depend heavily on specific hardware and drivers to be effective (even if you have the hardware and drivers, getting it to work on a VM is something different).

I've found a number of solutions in R, Python, and ImageJ which work well on the small scale or with specialized types of photos like medical images, satellite or microscopic images, but not for real world photos with lots of "noise" and features which can vary a good bit between photos.

Is there a Windows-friendly or platform-agnostic scalable solution that's free and minimally complicated?


问:
Pure math questions arising in computer vision, and the need for mathematicians to solve them
I've read in several articles/answer written by CS (grad) students or other people that the advanced knowledge of pure math that's required for certain areas of computer vision is sometimes too high for the CS student themselves, and requires the involvement actual math people.

I do wonder if it's true at all? I don't certainly pose this question from a negative view. As a new computer vision researcher who just switched from pure math, I find the math behind CV interesting, but the questions asked to be of different nature. But can the CV researchers with CS background not answer those math questions?


问:
What are alternatives to openbabel to convert CML file to image?
I recently started to perform natural resonance theory (NRT) analysis with NBO6. This program may automatically create a chemical markup language (cml) file with the treated structures. Unfortunately I don't speak cml well enough to judge if the file generated is any good.

What I would like to do is to convert this generated file into an image. That means, that the included 3D coordinates have to be flattened to 2D. I think the problem lies at this step as I get an erroneous structure reported.

As a test molecule I chose acetic acid and I am including a total of three resonance structures (It gets worse when I use more). The first and main contributor is displayed just fine. The second has missing formal charges. These are not explicitly given in the cml file, so they would need to be calculated. Openbabel is certainly able to do that, but for some reason, it is not working here. The third structure is completely ripped apart. The atoms should be at the same position as in the first and second image.

I am using a self-compiled local version of openbabel (2.3.2) with the following command (test.cml at the end):

```<code>  ~/local/tools/openbabel-install/bin/babel -icml test.cml -opng test.png
</code>```

I obtain the following image:

Along with the error messages:

```<code>==============================
*** Open Babel Error  in DoHCounts
  In atom a1 the number of explicit hydrogens exceeds the hydrogenCount attribute.
==============================
*** Open Babel Error  in DoHCounts
  In atom a1 the number of explicit hydrogens exceeds the hydrogenCount attribute.
==============================
*** Open Babel Error  in DoHCounts
  In atom a1 the number of explicit hydrogens exceeds the hydrogenCount attribute.
3 molecules converted
3 errors 109 audit log messages
</code>```

I don't have any idea how to access those audit log messages (/var/log/messages: Permission denied) and if they are in any way helpful.

However, the error messages already seem to indicate, that the generated cml file is faulty.

```<code><?xml version="1.0"?>
<cml title="Resonance forms from NRT analysis">
 <molecule id="1" title="Wgt=64.49%">
  <atomArray>
   <atom id="a1" elementType="C" hydrogenCount="0" x3="-1.40" y3="-0.11" z3="0.00" lonePair="0"/>
   <atom id="a2" elementType="H" hydrogenCount="0" x3="-1.93" y3="0.86" z3="0.00" lonePair="0"/>
   <atom id="a3" elementType="H" hydrogenCount="0" x3="-1.70" y3="-0.70" z3="0.89" lonePair="0"/>
   <atom id="a4" elementType="H" hydrogenCount="0" x3="-1.70" y3="-0.70" z3="-0.89" lonePair="0"/>
   <atom id="a5" elementType="C" hydrogenCount="0" x3="0.09" y3="0.13" z3="-0.00" lonePair="0"/>
   <atom id="a6" elementType="O" hydrogenCount="0" x3="0.78" y3="-1.05" z3="0.00" lonePair="2"/>
   <atom id="a7" elementType="O" hydrogenCount="0" x3="0.65" y3="1.21" z3="0.00" lonePair="2"/>
   <atom id="a8" elementType="H" hydrogenCount="0" x3="1.73" y3="-0.80" z3="0.00" lonePair="0"/>
  </atomArray>
  <bondArray>
   <bond atomRefs2="a1 a2" order="1"/>
   <bond atomRefs2="a1 a3" order="1"/>
   <bond atomRefs2="a1 a4" order="1"/>
   <bond atomRefs2="a1 a5" order="1"/>
   <bond atomRefs2="a5 a6" order="1"/>
   <bond atomRefs2="a5 a7" order="2"/>
   <bond atomRefs2="a6 a8" order="1"/>
  </bondArray>
 </molecule>
 <molecule id="2" title="Wgt=23.43%">
  <atomArray>
   <atom id="a1" elementType="C" hydrogenCount="0" x3="-1.40" y3="-0.11" z3="0.00" lonePair="0"/>
   <atom id="a2" elementType="H" hydrogenCount="0" x3="-1.93" y3="0.86" z3="0.00" lonePair="0"/>
   <atom id="a3" elementType="H" hydrogenCount="0" x3="-1.70" y3="-0.70" z3="0.89" lonePair="0"/>
   <atom id="a4" elementType="H" hydrogenCount="0" x3="-1.70" y3="-0.70" z3="-0.89" lonePair="0"/>
   <atom id="a5" elementType="C" hydrogenCount="0" x3="0.09" y3="0.13" z3="-0.00" lonePair="0"/>
   <atom id="a6" elementType="O" hydrogenCount="0" x3="0.78" y3="-1.05" z3="0.00" lonePair="1"/>
   <atom id="a7" elementType="O" hydrogenCount="0" x3="0.65" y3="1.21" z3="0.00" lonePair="3"/>
   <atom id="a8" elementType="H" hydrogenCount="0" x3="1.73" y3="-0.80" z3="0.00" lonePair="0"/>
  </atomArray>
  <bondArray>
   <bond atomRefs2="a1 a2" order="1"/>
   <bond atomRefs2="a1 a3" order="1"/>
   <bond atomRefs2="a1 a4" order="1"/>
   <bond atomRefs2="a1 a5" order="1"/>
   <bond atomRefs2="a5 a6" order="2"/>
   <bond atomRefs2="a5 a7" order="1"/>
   <bond atomRefs2="a6 a8" order="1"/>
  </bondArray>
 </molecule>
 <molecule id="3" title="Wgt=4.20%">
  <atomArray>
   <atom id="a1" elementType="C" hydrogenCount="0" x3="-1.40" y3="-0.11" z3="0.00" lonePair="0"/>
   <atom id="a2" elementType="H" hydrogenCount="0" x3="-1.93" y3="0.86" z3="0.00" lonePair="0"/>
   <atom id="a3" elementType="H" hydrogenCount="0" x3="-1.70" y3="-0.70" z3="0.89" lonePair="0"/>
   <atom id="a4" elementType="H" hydrogenCount="0" x3="-1.70" y3="-0.70" z3="-0.89" lonePair="0"/>
   <atom id="a5" elementType="C" hydrogenCount="0" x3="0.09" y3="0.13" z3="-0.00" lonePair="0"/>
   <atom id="a6" elementType="O" hydrogenCount="0" x3="0.78" y3="-1.05" z3="0.00" lonePair="3"/>
   <atom id="a7" elementType="O" hydrogenCount="0" x3="0.65" y3="1.21" z3="0.00" lonePair="1"/>
   <atom id="a8" elementType="H" hydrogenCount="0" x3="1.73" y3="-0.80" z3="0.00" lonePair="0"/>
  </atomArray>
  <bondArray>
   <bond atomRefs2="a1 a2" order="1"/>
   <bond atomRefs2="a1 a3" order="1"/>
   <bond atomRefs2="a1 a4" order="1"/>
   <bond atomRefs2="a1 a5" order="1"/>
   <bond atomRefs2="a5 a7" order="3"/>
   <bond atomRefs2="a6 a8" order="1"/>
  </bondArray>
 </molecule>
</cml>
</code>```

If you look at the coordinate information, you will see, that they are identical. The bond array just changes along with the lone pairs. It is also obvious from where the above mentioned errors stem from. Each entry has hydrogenCount="0", which is maybe correct because there are explicit hydrogens. But I am not sure. I found an open bug report referring to this attribute.

Now that I seem to not find a reason for this behaviour I need to find an alternative. The currently only viable and reliable method for me is to interpret the output of the NRT by hand.

Resulting from this I have a few questions:

Is the produced cml file correct? Does it contain errors and hence the problem is not with openbabel?

If the file is correct, are there alternatives to openbabel, which treat it correctly and which are able to convert it to a 2D image (structural representation)? (I would prefer a commandline tool.)

The output should be similar to this manually created image:

I have tried other programs, so far with not much luck.

Chemdraw can read the cml file and produces the structures with the correct bonds and the atoms are in the correct position. However, charges are missing and all molecules are at the same position, thus making it quite hard to obtain a good image. It also fails to read the title statement.

JChemPaint while technically able to read cml, fails to read this file. I could pinpoint it to the third structure, which produces an error and aborts. If only the first two are included it also fails to produce the charges.

ChemSketch cannot read cml formatting.


问:
Applying Newton-Raphson method to system of two differential equations, one time independent, one time dependent
The goal is to couple system of two nonlinear differential equations by applying appropriate space and time discretization and Newton-Raphson scheme.

The equations system is
$$\left[
\begin{array}{c}
F_{1}(u,v)       \\
F_{2}(u,v)
\end{array}\right]
=
\left[
\begin{aligned}
-\frac{\partial^2 u}{\partial x^2} &= v     \\
\frac{\partial v}{\partial t}
&=\frac{\partial J(u,v)}{\partial x}
\end{aligned}
\right]
$$
where $u(x,t), v(x,t)$ is function of time and space.

Notice $F_{1}$ is differential equation in space exclusively, and $F_{2}$ is partial differential equation in space and time.

The discretization in space and time are subscript $i$ and superscript $j$ respectively.

For $F_1$, discretizing in space gives
$$
F_{1}(u,v)\quad\rightarrow\quad
\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}+v_{i}=0
$$
The next is discretizing $F_2$ in space.
$$
F_{2}(u,v)\quad\rightarrow\quad
\frac{\partial v}{\partial t}
=
\frac{J_{i+1}-J_{i}}{\Delta x}
$$

The assumption is made; $J_{i+1}$ is function of $u_{i}$, $u_{i+1}$, $v_{i}$ and $v_{i+1}$ and  $J_{i}$ is function of $u_{i-1}$, $u_{i}$, $v_{i-1}$ and $v_{i}$.

In discretizing $F_{2}$ in time, the method TR-BDF2(Trapezoidal Rule - Second-order Backward Differentiation Formula) is applied, which is
$$
\frac{\partial q(z)}{\partial t}
=
f(z,t)
\quad\rightarrow\quad
\frac{q^{j+1}-q^{j}}{\Delta t}=
\frac{f^{j+1}+f^{j}}{2}
\quad\rightarrow\quad
2q^{j+1}-\Delta t f^{j+1}
-2q^{j}-\Delta t f^{j}
=0
$$

Once again, superscript $j$ and $j+1$ denotes discretization in time, not power.

Applying TR-BDF2 method to space discretized $F_{2}$ gives
$$
F_{2}(u,v)\quad\rightarrow\quad
2v_{i}^{j+1}- \frac{\Delta t}{\Delta x}(J_{i+1}^{j+1}-J_{i}^{j+1})
-2v_{i}^{j}- \frac{\Delta t}{\Delta x}(J_{i+1}^{j}-J_{i}^{j})
=0
$$

Now, there are two equations, one in only discretized in space and the other discretized in both time and space. How do these two equations be coupled in Newton-Raphson scheme?

Does applying TR-BDF2 to $F_{1}$ valid? Since $F_{1}$ has to remain invariant through time...

Anyway, $q$ is zero in $F_{1}$, so $f^{j+1}+ f^{j}=0$ (??) Then there are two equations both in time and space discretized.
$$
F_{1}(u,v)\quad\rightarrow\quad
\frac{u_{i+1}^{j+1}-2u_{i}^{j+1}+u_{i-1}^{j+1}}{\Delta x^2}+v_{i}^{j+1}
+\frac{u_{i+1}^{j+1}-2u_{i}^{j}+u_{i-1}^{j}}{\Delta x^2}+v_{i}^{j}
=0
$$
$$
F_{2}(u,v)\quad\rightarrow\quad
2v_{i}^{j+1}- \frac{\Delta t}{\Delta x}(J_{i+1}^{j+1}-J_{i}^{j+1})
-2v_{i}^{j}- \frac{\Delta t}{\Delta x}(J_{i+1}^{j}-J_{i}^{j})
=0
$$
Applying Newton-Raphson Scheme method to this system is,
$$
u^{k+1}=u^{k}+\delta u^{k},
\qquad\qquad
\left[
\mathbf{J}(u^{k})
\right]\delta u^{k}
=
-
\left[
\mathbf{F}(u^{k})
\right]
$$

where Jacobian and $\mathbf{F}$ is

$$
\left[
\begin{array}{*{8}c}
U_{1D} & U_{1U} &        &       & V_{1D} &      &      &       \\
U_{1L} & U_{1D} & U_{1U} &       &        &V_{1D}&      &       \\
       & U_{1L} & \ddots & U_{1U}&        &      &\ddots&      \\
       &        & U_{1L} & U_{1D}&        &      &      &V_{1D}\\
U_{2D} & U_{2U} &        &       & V_{2D} &V_{2U}&      &         \\
U_{2L} & U_{2D} & U_{2U} &       & V_{2L} &V_{2D}&V_{2U}&          \\
       & U_{2L} & \ddots & U_{2U}&        &V_{2L}&\ddots&V_{2U}  \\
       &        & U_{2L} & U_{2D}&        &      &V_{2L}&V_{2D}\\
\end{array}
\right]
\left[
\begin{array}{c}
\delta u_2\\
\vdots \\
\delta u_{n-1}\\
\delta v_2\\
\vdots \\
\delta v_{n-1}
\end{array}
\right]
=
-
\left[
\begin{array}{c}
F_1(i=2)\\
\vdots\\
F_1(i=n-1)\\
F_2(i=2)\\
\vdots\\
F_2(i=n-1)
\end{array}
\right]
$$

where $F_{1,2}(i)$ vector is inserting current $i$ value for $j$ terms and guess value $i$ for $j+1$ (next time) terms.

and each components are (D=main diagonal, U=first upper diagonal, L=first lower diagonal)

$$
\begin{aligned}
U_{1D}
&=
\frac{d}{du_i^{j+1}}F_{1}^{j+1}
=-\frac{2}{\Delta x^2}
& (2\leq i\leq n-1)
\\
U_{1U}
&=
\frac{d}{du_{i+1}^{j+1}}F_{1}^{j+1}
=\frac{1}{\Delta x^2}
& (2\leq i\leq n-2)
\\
U_{1L}
&=
\frac{d}{du_{i-1}^{j+1}}F_{1}^{j+1}
=\frac{1}{\Delta x^2}
& (3\leq i\leq n-1)
\\
V_{1D}
&=
\frac{d}{dv_{i}^{j+1}}F_{1}^{j+1}
=1
& (2\leq i\leq n-1)
\\
U_{2D}
&=
\frac{d}{du_i^{j+1}}F_{2}^{j+1}
=-\frac{\Delta t}{\Delta x^2}\left(
\frac{dj_{i+1}^{j+1}}{du_{i}^{j+1}}
-
\frac{dj_{i}^{j+1}}{du_{i}^{j+1}}
\right)
& (2\leq i\leq n-1)
\\
U_{2U}
&=
\frac{d}{du_{i+1}^{j+1}}F_{2}^{j+1}
=-\frac{\Delta t}{\Delta x^2}\left(
\frac{dj_{i+1}^{j+1}}{du_{i+1}^{j+1}}
\right)
& (2\leq i\leq n-2)
\\
U_{2L}
&=
\frac{d}{du_{i-1}^{j+1}}F_{2}^{j+1}
=\frac{\Delta t}{\Delta x^2}\left(
\frac{dj_{i-1}^{j+1}}{du_{i-1}^{j+1}}
\right)
& (3\leq i\leq n-1)
\\
V_{2D}
&=
\frac{d}{dv_i^{j+1}}F_{2}^{j+1}
=2-\frac{\Delta t}{\Delta x^2}\left(
\frac{dj_{i+1}^{j+1}}{dv_{i}^{j+1}}
-
\frac{dj_{i}^{j+1}}{dv_{i}^{j+1}}
\right)
& (2\leq i\leq n-1)
\\
V_{2U}
&=
\frac{d}{dv_{i+1}^{j+1}}F_{2}^{j+1}
=-\frac{\Delta t}{\Delta x^2}\left(
\frac{dj_{i+1}^{j+1}}{dv_{i+1}^{j+1}}
\right)
& (2\leq i\leq n-2)
\\
V_{2L}
&=
\frac{d}{dv_{i-1}^{j+1}}F_{2}^{j+1}
=\frac{\Delta t}{\Delta x^2}\left(
\frac{dj_{i-1}^{j+1}}{dv_{i-1}^{j+1}}
\right)
& (3\leq i\leq n-1)
\end{aligned}
$$

Again,  $J_{i+1}$ is function of $u_{i}$, $u_{i+1}$, $v_{i}$ and $v_{i+1}$ and  $J_{i}$ is function of $u_{i-1}$, $u_{i}$, $v_{i-1}$ and $v_{i}$.

Is this correct way?


问:
Need help writing the code for the following optimization
I need to find $X$ for which the following expression is minimized:

$$ \min ||Y - F^{-1}(X)||_2 + ||X||_1 $$

where $Y$ is an array (of length approx 44000, an audio sample I will be reading using Matlab's audioread( ) function) and $F^{-1}$ is the inverse fast Fourier transform ifft( ). Subscripts indicate 2-norm and 1-norm respectively.
Could someone help me going about implementing it using Matlab?


问:
Unwanted Oscillation in FDM simulation of elastic wave equation
I am using staggered grid FDTD for solving elastic wave equation. A description of which can be found at (geodynamics.usc.edu/~becker/teaching/557/reading/Virieux1987.pdf). I have generated a synthetic seismogram (description below) using it.

Synthetic seismogram : To generate a synthetic seismogram, a source is placed at top and an array of receiver is laid along the surface. A synthetic seismogram can be visualize as in this youtube animation (www.youtube.com/watch?v=_Gjr5EWNIQk)

My Problem is the spurious oscillations generated in FDTD simulation for elastic wave. I have taken care of grid dispersion which requires 10 grid points per wavelength (minimum). still I feel there is a lot of dispersion in it.
The model I have taken is a homogeneous elastic model with absorbing (damp) type  boundary conditions along all the edges.

Vp = 2000 m/s; Vs = 1500 m/s;

Grid spacing: dx = dv = 5m;

time sampling: dt = 1msec

Source: Ricker (20Hz).

Following is the response of at one receiver (k/a trace) and a snap of wave

EDIT 1

Presently I don't have privileged to add more than two link so I edited and joined two pics of traces for soft and hard types of sources.

The corrected trace after applying pressure source (soft type) is shown now. But still there are wriggles trailing the wave


问:
Would anyone give some review on adjoint methods in shape optimization method <in fluid dynamics>?

Modification

Is a continuous adjoint approach a reliable/promising method for optimizing the shape of an airfoil under conditions of unsteady turbulence flow?

Original Question

I saw lots of people is doing shape optimization of airfoil using optimization theory.

But my feeling is that its effectiveness is very limited, and the formulation in continuous adjoint and programming in discrete adjoint method are all very complicated.

And I do not think it is very promising after a few month's tests and programming. It seems the "improved" shape of surface is very noisy. And many of its assumptions only valid in frozen turbulence.

Could anyone share more his experience on this subject?

Thanks

OwnerDisplayName: Daniel


问:
Does scaling factor affect discretization?
Suppose I want to solve the below equation numerically.

$$
\frac{dy}{dx}=y
$$

I'd like to normalize the space discretization by choosing
$$
a\bar{x}=x
$$
where I assume $\bar{x}$ is unity.
Then the equation is

$$
\frac{1}{a}\frac{dy}{d\bar{x}}=y
\qquad\rightarrow\qquad
\frac{dy}{d\bar{x}}=ay
$$

Thus, my discretization becomes,

1)

$$
\frac{y_{i+1}-y_{i}}{\Delta \bar{x}}=ay
$$

which I don't think it makes sense, since $\Delta \bar{x}$ must be unity, so is

2)

$$
y_{i+1}-y_{i}=ay
$$

right?


问:
How to solve bring my implicit equation to closed form?
I have a simulated system as shown in the following image:

$L_0$ is attached to two other bodies $L_1$ and $L_2$. Furthermore, body $L_3$ is also in the simulation (it is attached to $L_2$ but that is not important). Each body has a position, so $L_0$'s position is $\vec{P}^0 = (x_0, y_0)$, $L_1$'s position is $\vec{P}^1 = (x_1, y_1)$, and so on.

$L_0$'s movement is governed by two forces, a linear attraction force to bodies it is attached to:

$\vec{A}_0 = \frac{\Delta P}{\Delta t} = k(\vec{P}^1 - \vec{P}^0) + k(\vec{P}^2 - \vec{P}^0)$

which adheres to Hook's Law of springs: $\vec{F} = -kr$ where $r$ is the distance between the two bodies. and a non-linear repulsion force which acts on $L_0$ from all other bodies:

$\vec{R} = \frac{\Delta P}{\Delta t} = \vec{U}^{01} \frac{Gm^0m^1}{|\vec{P}^0 - \vec{P}^1|^2} + \vec{U}^{02} \frac{Gm^0m^2}{|\vec{P}^0 - \vec{P}^2|^2} + \vec{U}^{03} \frac{Gm^0m^3}{|\vec{P}^0 - \vec{P}^3|^2}$

Where $\vec{U}^{01}$ is the direction vector from $L_1$ to $L_0$, so

$\vec{U}^{01} = \frac{\vec{P}^1 - \vec{P}^0}{|\vec{P}^1 - \vec{P}^0|}$ and the same is true for $\vec{U}^{02}$ and $\vec{U}^{03}$

This force adheres to Newton's law of gravity between two bodies with masses $m^1$ and $m^2$, $\vec{F} = \frac{Gm^1m^2}{r^2}$ where $r$ is the distance between the two bodies (so the repulsion increases quadratically as the bodies get closer)

In essence then, $L_0$'s movement is governed by the equation:

$\frac{\Delta P}{\Delta t} = k(\vec{P}^1 - \vec{P}^0) + k(\vec{P}^2 - \vec{P}^0) + \vec{U}^{01} \frac{Gm^0m^1}{|\vec{P}^0 - \vec{P}^1|^2} + \vec{U}^{02} \frac{Gm^0m^2}{|\vec{P}^0 - \vec{P}^2|^2} + \vec{U}^{03} \frac{Gm^0m^3}{|\vec{P}^0 - \vec{P}^3|^2}$

In my simulation, I want to use the implicit (Backward Euler) method of the discretized equation, so I discretize to:

$\vec{P}_{n+1}^0 = \vec{P}_n^0 + \Delta tk(\vec{P}_n^1 - \vec{P}_n^0) + \Delta tk(\vec{P}_n^2 - \vec{P}_n^0) + \Delta t\vec{U}_n^{01} \frac{Gm^0m^1}{|\vec{P}_n^0 - \vec{P}_n^1|^2} + \Delta t\vec{U}_n^{02} \frac{Gm^0m^2}{|\vec{P}_n^0 - \vec{P}_n^2|^2} + \Delta t\vec{U}_n^{03} \frac{Gm^0m^3}{|\vec{P}_n^0 - \vec{P}_n^3|^2}$

and then substitute $P_{n+1}^0$ instead of $P_n^0$ at the appropriate places to get:

$\vec{P}_{n+1}^0 = \vec{P}_n^0 + \Delta tk(\vec{P}_n^1 - \vec{P}_{n+1}^0) + \Delta tk(\vec{P}_n^2 - \vec{P}_{n+1}^0) + \Delta t\frac{\vec{P}_{n+1}^0 - \vec{P}_n^1}{|\vec{P}_{n+1}^0 - \vec{P}_n^1|} \frac{Gm^0m^1}{|\vec{P}_{n+1}^0 - \vec{P}_n^1|^2} + \Delta t\frac{\vec{P}_{n+1}^0 - \vec{P}_n^2}{|\vec{P}_{n+1}^0 - \vec{P}_n^2|} \frac{Gm^0m^2}{|\vec{P}_{n+1}^0 - \vec{P}_n^2|^2} + \Delta t\frac{\vec{P}_{n+1}^0 - \vec{P}_n^3}{|\vec{P}_{n+1}^0 - \vec{P}_n^3|} \frac{Gm^0m^3}{|\vec{P}_{n+1}^0 - \vec{P}_n^3|^2}$

This is where I'm stuck. How do I solve for $P_{n+1}^0$? Can I isolate it? do I need to use Newton's method, and if so how do I do that for a 2-coordinate vector?

The example talks about $L_0$, but I actually have to do it for each of the bodies described. So I hope the answer I get could be generalized for $L_1$, $L_2$, and $L_3$


问:
Quasi Newton method for block diagonal Hessian
I am having an unconstrained optimization problem, where the Hessian is positive semidefinite and block diagonal. The function is strictly convex, hence, the curvature condition ($ s^{T}_{k}y_{k} > 0$) is supposed to be always satisfied. Here, $s_{k}$ and $y_{k}$ are the usual secant pairs, where $s_{k} = x_{k}-x_{k-1}$ and $y_{k}=g_{k}-g_{k-1}$, where $g_{k}$ is the gradient at iteration $k$.

I am noticing that this curvature condition is indeed getting satisfied if I take the entire variable vector and gradient vector. But I was hoping that it will be satisfied at the block level, also, for all the blocks. i.e., by restricting $s$ and $y$ to block level vectors, for all the blocks. This is not being observed for all the blocks. My intuition says that it should happen because the blocks of a PSD matrix are also PSD. Can somebody give more understanding about this?


问:
Concave Cells and Finite Volume Method
I would like to know how to compute the "best" cell center for pressure calculation in case of concave polyhedral cells, if it makes any sense. My face centers are defined as faces barycenters and I'm just wondering how to correctly determine the pressure-point.

Any hint?

Title: Concave Cells and Finite Volume Method


问:
I'm using linear programming for production planning. Does the order in which I make products affect the cost?
I have a collection of different scrap aluminium alloys. I want to mix them together to make new alloys with customer-defined compositions.

Sometimes this will involve little more than melting down several tonnes of scrap aluminium cans to make ingots, which are then sold to manufacturers of aluminium cans so that they can make more aluminium cans.

Sometimes the required product will be very different to any of the relatively cheap scrap metal. This will require mixing different scrap metals with relatively expensive additions (known in the trade as "master alloys").

Here's the problem.

I schedule production in the order A, B, C, D and apply the simplex algorithm to each product using Excel solver. The first time I do this, there will be a certain number of tonnes of each type of scrap. The algorithm will calculate the cheapest mixture of scrap and master alloys to make product A.

However, when I do the same for product B, I have to deduct the quantities of scrap and master alloys used whilst making product A (I can't use the same material twice!). I clearly have to do the same when I make product C, then product D.

In other words, my material constraints change each time I make a product. And the sequence of material constraints will be different for each production sequence. The sequence A, B,C, D will give a different order of material constraints to the order D, C, B, A.

I assume that I can find the cheapest total cost simply by aggregating the whole day's production into one huge block of metal and running Excel solver on that. However, that doesn't tell me how to achieve that minimum cost in a real production schedule.

I have N products to make in any given day. Is there a way of finding the globally cheapest production schedule order without trying out all N! permuatations?


问:
Minimization of the sum of convex function and non-linear non convex function
I'm trying to minimize the unconstrained scalar sum of a quadratic convex function (to which a convex optimizer is readily applied) and a non-linear and non-convex function which is differentiable.

General purpose non-linear optimizers are not providing adequate results, which I assume is due to the high dimensionality of my problem (>=10000), so I was wondering if there was a better way to go about minimizing my function given that it is decomposable into a convex and differentiable non convex part.


问:
Fixing a near singular covariance matrix
Given a near singular covariance matrix, the standard method of 'fixing' it seems to be to add a small damping coefficient $c>0$ to the diagonal, which serves to bump all the eigenvalues up by this amount.

Is there a way to determine the optimal value of this damping coefficient without resorting to some iterative procedure where I jump back and forth between inverting the matrix to look for infinities and NaNs and making $c$ increasingly small?

I could look at the minimum eigenvalue fairly cheaply, but it seems like the correct choice for $c$ would still also depend on both the entries of the matrix as well as its size.

Finally, it's more important that my algorithm always works than it is that the perturbation is as small as possible, so I'm willing to accept more error for the sake of a procedure with some sort of guarantee that it will always work.


问:
Finding errors in frequency from a Fast Fourier Transform from Gaussian fitting
I took a FFT of sound in a box generated by a frequency sweep over a range of frequencies, and have an array of frequencies and their corresponding FFT amplitudes.  According to models for the acoustic waves, certain frequencies should be resonant frequencies, which certainly shows in the data as Gaussian-like peaks.

I was wondering if there might be some program or algorithm that can quickly determine the peaks in the spectrum as well as their standard deviations (if assumed to be Gaussian).

I already have a brute force way of doing, that is to check each point and their values on Igor Pro, and finding the value of $f$ such that the spectrum amplitude, $A = A_0 e^{-1/2}$ where $A_0$ is the peak.

This is due to the fact that the equation of the Gaussian is
$$
A_0 e^{-\frac{f^2}{2\sigma_f}}
$$
where $\sigma_f$ is the standard deviation in $f$, so that at $f=\sigma_f$, the exponential is $e^{-1/2}$.

I was wondering if there was a way to automate this.  Finding the peaks themselves would probably involve some sort of discrete derivative to find the extrema (which would have to ignore turning points due to noise fluctuations).  Then perhaps I could compare values over a range about the peaks, but of course the comparisons won't be exact equality, but instead finding the closest values, since I am working with discrete data.


问:
Solving Schrodinger's Equation Numerically in a Bunimovich Stadium
I need to solve, as mentioned, Schrodinger's equation in a Bunimovich stadium-shaped infinite potential well with Dirichlet BC Numerically (this isn't possible analytically). In order to do so, I need to define the discrete Hamiltonian matrix of the system. Since the shape isn't defined well in any coordinates system, I have no idea how to do so.

Normally (as I did with a square well), I'd define the size of the well and then solve 1D TISE for each dimension using MATLAB eig function. Can it be solved similarly?

Any help will be greatly appreciated.


问:
Beowulf Cluster For MD simulations
What would you recommend as an all-in-one cluster management software for a Beowulf cluster? We have about 7(most of them are 8 core machines and 1 is a 48 core machine) very good compute nodes in our research group and they are currently managed using Rocks and run Centos 6. We have been doing MPI-based computations, particularly, molecular simulations (Molecular Dynamics and Monte Carlo Simulations) reliably but we would now like to update the software. I suppose doing a system upgrade would break packages and reinstalling the cluster is a better option. I am wondering if we should use Rocks again or go for some other all-in-one solutions. We would want something that's easy to maintain (we all have experience with linux but not server management) and has LTS.

So, my questions are: 1. Is there a Debian-based option like PelicanHPC that is good for deploying on a long term basis? 2. Is building a beowulf cluster based on Debian/Ubuntu/Centos 7 from scratch a good idea?

Title: Beowulf Cluster For MD simulations


问:
Very long running time for Haase and Muller (2014) coded in Python + Gurobi
I have read the paper of Haase and Muller (2014) where they present the linear reformulations of the multinomial logit choice probabilities and compare different approaches. I have tried for an exercise to code their models, particularly the variant of Haase (2009), and to solve it with Gurobi 6.5. The issue is that I obtain really long running times.

In their numerical investigation they used CPLEX 12 on a 64-bit Windows Server 2008 with 1 Intel Xeon 2.4 GHz processor and 24 GiB of RAM. On the other hand I am using Python + Gurobi on a 64-bit Windows 8 with 2 Six-Core AMD Opteron 2.00 GHz processors and 32 GiB of RAM.

I understand that there can be differences between Gurobi and CPLEX in running times on different models/programs. In some cases the Gurobi is better, and in some cases the CPLEX is better. Also, I understand that my machine is not very representative and I hope that 0.4 GHz doesn't make such a big difference. But the difference is really huge. For some small problems, that Haase and Muller needed a little bit more than a minute to solve them, on my machine I need more than two hours. Is it possible that choice of the solver and machine can matter that much? Can someone, please, check my code, because I am newbie in this game.

This is the mathematical model written in LaTeX:

$ \max F = \sum_{i\in I} \sum_{j\in J} x_{ij} $

s.t.

$ \sum_{j\in J} y_j = r $

$ \tilde{x}_{i} + \sum_{j\in J} \leq 1, \forall i \in I $

$ x_{ij} - \frac{\phi_{ij}}{1 + \phi_{ij}} y_j \leq 0, \forall i \in I, j \in J $

$ x_{ij} - \phi_{ij}\tilde{x}_i \leq 0, \forall i \in I, j \in J $

$ y_i \in \{0, 1\}, \forall j \in J $

$ x_{ij} \geq 0, \forall i \in I, j \in J $

$ \tilde{x}_i \geq 0, \forall i \in I $

where $ \phi_{ij} = \frac{ e^{v_{ij}} }{\sum_{k\in M\setminus J} e^{v_{ik}}} $.

Here is the code in Python.

```<code>from gurobipy import *
import math
import numpy as np

n = 200 # The number of customers
m = 60 # The number of possible facility locations.
x_coords = 30 * np.random.random_sample(n + m)
y_coords = 30 * np.random.random_sample(n + m)
I = { i: (float(x_coords[i]), float(y_coords[i])) for i in range(n) } # The customers' locations.
M = { j: (float(x_coords[n + j]), float(y_coords[n + j])) for j in range(m) } # The possible facility locations.
J = {i: M[i] for i in range(m-10)} # J is subset of M. The competitor has located facilities in M\J.
MdiffJ = {i: M[m-10+i] for i in range(10)} # Competitor's facility locations. |M\J| = 10 as in the paper.
d = {}
v = {}
for i in I:
    for j in J:
        d[i,j] = abs(I[i][0] - J[j][0]) + abs(I[i][1] - J[j][1]) # Rectangular distances.
        v[i,j] = -0.2 * d[i,j]
r = round(0.2 * len(J)) # == 0.2 * |J|, as in the paper.

model = Model()

# Creating variables
y = {}
x_tilda = {}
x = {}
for j in J:
    y[j] = model.addVar(vtype=GRB.BINARY, name="y_%i" % j)
for i in I:
    x_tilda[i] = model.addVar(lb=0, name="x_tilda_%i" % i)
    for j in J:
        x[i,j] = model.addVar(lb=0, name="x_%i%i" % (i,j))

# Integrating variables
model.update()

# Setting the objective.
obj = LinExpr()
obj += quicksum(quicksum(x[i,j] for j in J) for i in I)

# Setting the objectiv as a maximization problem.
model.setObjective(obj, GRB.MAXIMIZE)

# Computing phi (constants)
phi = {}
for i in I:
    for j in J:
        phi[i,j] = math.e**v[i,j] / (sum(math.e**v[i,k] for k in MdiffJ))

# Adding constraints.
model.addConstr(quicksum(y[j] for j in J) - r == 0, "sum(y_*) == r")
for i in I:
    model.addConstr(x_tilda[i] + quicksum(x[i,j] for j in J) <= 1, "x_tilda_%i+sum(x_%i*) <= 1" % (i,i))
    for j in J:
        model.addConstr(x[i,j] - (phi[i,j] / (1 + phi[i,j])) * y[j] <= 0, "x_%i%i - frac * y_%i <= 0" % (i,j,j))
        model.addConstr(x[i,j] - phi[i,j] * x_tilda[i] <= 0, "x_%i%i - phi_%i%i * x_tilda_%i <= 0" % (i,j,i,j,i))

# Computing the optimal solution.
model.optimize()
</code>```


问:
Initial conditions of nbody problem
Note: I have edited the original question to try to be clearer.

I'm trying to simulate a gravitational nbody problem with a Barnes-Hut algorithm, in which the bodies start in a homogeneous sphere with random velocities. Because of this, the initial angular momentum is almost zero. I have tested the code with things like the solar system and compared its results to a brute force algorithm and it seems to work fine.

Despite using very low time steps (I'm using a velocity verlet integrator, if that matters) the angular momentum increases a lot, not being conserved at all. This image shows a typical increase of the momentum (later I will show more detailled data)

You can see that the energy is being conserved perfectly but the angular momentum keeps increasing with the bodies being accelerated. All the angular momentums are computed from the actual center of mass of the cluster.

When a body has a close encounter and gets expelled from the cluster, both the angular momentum and the energy increase a lot:

I have tried different timesteps (from $\delta t= 10^{-8}t_0$ to $\delta t=10^{-3}t_0$) and gravitational softening parameters $\epsilon$ (from $F_{ij} \propto r_{ij}(\epsilon + r_{ij}^3)^{-1}$) from $1$ light year to $10^{-5}$ light years (the initial separation between bodies is $6$ ly), and they all have the same general results (there are 8 series in each line but the deviation is minimal):

The sudden rise of the lines is caused by a body being expelled by a close encounter. My principal concern is the initial increase of the angular momentum, because it doesn't come with a variation in the total energy. Because real galaxies have angular momentum, I thought that maybe I should start the simulation with a non zero angular momentum.

My questions are:

What are the standard initial conditions in this kind of problems? They involve a small initial angular momentum or they use a "spinning" cluster to begin the simulation?

Which variation in angular momentum and energy can be considered reasonable in nbody simulations?

Edit: I wasn't aplying gravitational softening, but it does not change the results.

Edit: This is the relevant code of the angular momentum, I can see a flaw but who knows...

```<code>double acc = 0;
for(int i=0;i<nstars;i++){
    double r[3] = { x[i]-cmx , y[i]-cmy , z[i]-cmz };
    double v[3] = { vx[i] , vy[i] , vz[i] };
    double li[3] = { r[1]*v[2]-r[2]*v[1], -(r[0]*v[2] - r[2]*v[0]),  r[0]*v[1] - r[1]*v[0] };
    // Li = Rᵢ×Vᵢ
    acc += m[i]*norm(li);
    // acc += m⋅R×V
    }
return acc;
</code>```

cmx, cmy, cmz are the coordinates of the mass (just checked, it doesn't move a bit and it's where it should be). Using a fixed reference frame makes no difference.

Edit: Just tried gravitational softening (it was badly implemented). It helps a lot to avoid bodies fly away, but makes no difference in the constant increase of the angular momentum.

Title: Initial conditions of nbody problem


问:
Phase dislocations and numerical accuracy
I am solving the nonlinear Schrodinger equation (NLSE),
$$A_t+iA_{xx}+i|A|^2A=0$$
where $A$ is a complex valued function, which can be written as $A=ae^{i\theta}$ for $a,\theta$ real. Now, for particular sets of initial conditions, I find that phase dislocations occur, and hence there are singularities in the frequency, $\omega=-\theta_t$.

My question is, do phase singularities effect the accuracy of the model? Furthermore, physically, why can rapid variations in the phase occur in an equation that is derived under the WKB assumption? Does this mean the equation is not uniformly valid?

An example will be provided below.

I am using a split-step pseudo spectral method. My code is written in Fortran but I have found I can reproduce the results with a Matlab code I found online, which is provided below.

The initial condition takes the form
$$A=sech(x)e^{iCx^2}$$

for $C$ a constant representing a chirp to the signal, which can lead to packet focusing due to dispersion. The evolution of the packet is shown in figure 1.

Now, the phenomena in question occurs around $t=0.20$ when the amplitude $a$ goes to 0, and hence the phase is not well defined. This leads to singularities in $\omega$. A better picture of this is shown in figure 2.

We see that as the amplitude goes towards 0 at the front and rear of the group, i.e. the packet splits, the phase starts to become singular. Note, I've decimated $\omega$ by a factor 10 to make my point.

Note, the amplitude $a$ heuristically seems to vary slowly, but the phase clearly has rapid variations where the packet amplitude becomes small. Does this mean my results are no longer physically? Finally, and this is off topic, what sets the sign of $\omega$ at these singularities?

Note, to check the numerical accuracy of the scheme I confirm the energy is conserved to one part in $10^6$. Furthermore, there exist exact solutions to the NLSE, which I confirm stay solutions for long times with the resolution and step size used. Finally, the CFL condition for linear wave speeds ($u=1$ here) is obeyed. I have played around with the resolution and time steps, and find they give the same results.

Matlab script:

```<code>  %  This Matlab script file solves the nonlinear Schrodinger equation
  %  with the first order split-step Fourier method due to Hardin and Tappert.
  %                              Andre Weideman, 1995
  %
  N=512; %number of fourier modes
  dt = 0.001;  %time step
  M=300; % number of time steps
  t=0:dt:M*dt-dt;
  L=4*pi; % size of the domain
  I=sqrt(-1);
  % J=100;
  h = L/N;                       % Space step
  n = [-N/2:1:N/2-1]';           % Indices
  x = n*h;            % Grid points
  %----------------------------------------
  % IC: Chirped sech packet
  u=sqrt(2)*sech(x).*exp(-I*x.^2);
  %----------------------------------------
  U = u;           % Compute initial condition; save it in U
  e = -4*n.*n*pi*pi/L/L;         % Squares of wavenumbers.
  for m = 1:1:M                  % Start time evolution
  u = exp(dt*I*(abs(u).*abs(u))).*u;  % Solve nonlinear part of NLS
  c = fftshift(fft(u));                 % Take Fourier transform
  c = exp(dt*I*e).*c;                   % Advance in Fourier space
  u = ifft(fftshift(c));                % Return to physical space
  U = [U u];
  end
</code>```


问:
Is casadi suitable for data fitting?
Quite often I do fit some ODE or DAE systems to my data (small to medium sized problems). Via the assimulo package, I found Casadi and read a bit about the language modellica.

Casadi offers automatic differentiation and an interface to ipopt. As far as I understood. It is designed to solve and optimize ODE or DAE systems but until now I havent found a report of someone using it for datafitting.

Would it be reasonable to use casadi (or jmodellica) for such a task ? The combination of Ipopt with Casadi seems to be very powerfull.

Title: Is casadi suitable for data fitting?


问:
Using entropy functions for increasing numerical stability
Regarding the numerical stabilization of two-dimensional advection equation,
\begin{equation}
\dfrac{\partial f}{\partial t} + \Big(\dfrac{d\varepsilon_1(k)}{dk}\Big)\dfrac{\partial f}{\partial z}
- \Big( \dfrac{d\varepsilon_2(z)}{dz} \Big)\dfrac{\partial f}{\partial k} = 0
\end{equation}
I have encountered some texts that propose using entropy functions. For example, it is mentioned that the choice of
\begin{equation}
h(f,z,k) = fe^{\varepsilon_1+\varepsilon_2}
\end{equation}
will do the job for the above equation.

How is the given function (or any other entropy function) increasing the numerical stability? What is the main idea behind using such functions?


问:
Optimization on the manifold of stochastic matrices
So I have an optimization problem of the form
$$\text{maximize}\hspace{3mm}f(A):{\bf R}^{K\times K}\rightarrow{\bf R}$$
$$\text{subject to}\hspace{19mm}A^T{\bf 1}=\bf{1}$$
$$\hspace{33mm}A\geq 0$$
where the constraint set is thus the manifold of left stochastic matrices (each column sums to 1).  The objective function $f$ is not concave and thus my goal is to very quickly find a local maxima (speed is important).

The function $f$ is a nasty piece of work, nevertheless the gradient $\nabla f$ is analytically computable.  However even without the constraint, setting the gradient equal to zero and solving for the extrema is intractable, and thus I can't imagine that back solving for Lagrange multipliers is doable.

The method I've tried is manifold gradient ascent using some code found here from the Manopt package.  They refer to the manifold of left stochastic matrices as the multinomial manifold, and the code provides methods for basically projecting the gradient onto the tangent space, and then retracting this new gradient onto the manifold itself.

Unfortunately as far as retractions go it's fairly primitive I think, and it involves entry-wise dividing the projected gradient by the current point and then exponentiating it, which causes overflow unless you make your step-sizes quite small (however I don't fully understand the intricacies of manifold optimization so I could be wrong about this).

I could compute the Hessian numerically and try to run a 2nd order optimization on the multinomial manifold, but I'd have to do some reading to figure out exactly how to do this.  Basically I need something that is fast, speed is more important than accuracy, as I have to do this thousands of times and it's only 1 step in a larger coordinate ascent algorithm.

The gradient $\nabla f$ is complicated and thus each evaluation of it is quite costly, however the one saving grace is that the dimension of $A$ is relatively small, almost certainly $<100$.

What would be the best approach here?

Update

As requested, here is the function and its gradient:
\begin{align}
f(A )&= \Big[-\psi\big(\sum_{i=1}^KC_{ij}\big)+\sum_{i=1}^KA_{ij}(\psi(C_{ij})-\log C_{ij})\Big]_{j=1:K} \cdot \big(\sum_{n=2}^NA^{n-2}\big)\cdot {\bf x}\\
&+ \sum_{n=1}^N[\log B_{in}]_{i=1:K}^T\cdot A^{n-1}\cdot{\bf x}
\end{align}

\begin{align}
\nabla f(A) &= [\psi(C_{ij}) - \log A_{ij} -1]_{i=1:K,\hspace{1mm}j=1:K}\cdot \text{diagm}\Big(\big(\sum_{n=2}^NA^{n-2}\big)\cdot{\bf x}\Big)\\
&+ \sum_{n=3}^{N}\sum_{r=0}^{n-3}\Big((A^r)^T\cdot\Big[-\psi\big(\sum_{i=1}^KC_{ij}\big)+\sum_{i=1}^KA_{ij}(\psi(C_{ij})-\log C_{ij})\Big]_{j=1:K}^T \cdot {\bf x}^T\cdot (A^{n-3-r})^T\Big)\\
&+ \sum_{n=2}^N\sum_{r=0}^{n-2}(A^r)^T\cdot[\log B_{in}]_{i=1:K}\cdot{\bf x}^T\cdot(A^{n-2-r})^T
\end{align}
Where $\psi$ is the digamma function, $A$ and $C$ are $K\times K$ matrices, $B$ a $K\times N$ matrix and ${\bf x}$ is a length $K$ vector, all real-valued.  The function $\text{diagm}$ converts a vector to a diagonal matrix.  Also I checked my derivation of the gradient numerically, so you can be sure its correct.

I'm writing this in Julia, and I tried using the NLopt package both using a gradient-based method and a derivative free method, but they were even slower than my original manifold approach.


问:
Coupled system of differential equations with eigenvalue
I have the following system of coupled differential equations (the Bogoliubov-de Gennes equations for a certain geometry):

$$
\begin{cases}
-\frac{1}{2}u''(r) + \big[\frac{1}{2r^2} + f_1(r)\big]u(r) + f_2(r)v(r)=\omega u(r)\\
-\frac{1}{2}v''(r) + \big[\frac{1}{2r^2} + f_1(r)\big]v(r) + f_2(r)u(r)=-\omega v(r),
\end{cases}
$$
where $f_i(r)$ are known functions which go to zero for large $r$ and are finite for small $r$. This is in fact an eigenvalue problem, because I need to find the value of $\omega$ for which this system has (normalizable) solution.

Now I would like to solve this system numerically: hence find the value of $\omega$ such that the system has a solution and then also the functions $u(r),v(r)$. I am familiar with numerical recipes to solve coupled differential equations, but now that this is an eigenvalue problem as well, I am at a loss. Could anyone give me some pointers on how to tackle such a problem? Some references where this is described would be appreciated as well!


问:
Compute sparsity pattern of $A^2$
Suppose we have a sparse matrix $A$.  Is there any way to compute just the sparsity pattern of $A^2 = A \cdot A$ (I do not actually need to know what exactly the nonzero value are) faster than to compute the full sparse product $A \cdot A$? </p>

Title: Compute sparsity pattern of $A^2$


问:
Tucker factorisation to compare multiple PCA decompositions?
This is an entry-level question for multiway matrix decompositions. I have a set/population $k$ of entities (here biological cells) for each of which I also have a number ($l$) of flavours of length $N$ vectors (genetic-like measurements for which $N \gg k$) which characterise them. Now if I take just one type of vector, and look at its values across all of the entities (cells), I do a PCA since I effectively have a $k \times N$ (depending on convention) array which I can explore in the conventional. If I repeat this analysis $l$ times for each flavour of the vector, I get $l$ very similar analyses but I would like to take an overall view. My question is this:

is it possible (or reasonable) to approach this as a multiway analysis on an $l \times k \times N$ data source says using a Tucker factorisation/tensor decomposition?

is this the most appropriate approach?

are there any available resources, papers or books that I should read to understand this?

The aim is to understand how the different vector flavours (types of genetic-like information) differ in the characterisation of the cells.

I'm currently working with the rTensor package in R which I think may be based upon a similar one in MATLAB. I have a maths background and am in principle not afraid of tensors although this would be the first time I have used them for data analysis.


问:
Solving a large system of nonlinear equations, where timeseries are the unknown
I am trying to solve a problem, which I find quite hard, like, headache-hard. I have to solve the following set of $M$ nonlinear equations:
$$F(X)=\begin{bmatrix}f_1 (X)\\f_2 (X)\\...\\f_M (X)\\ \end{bmatrix}=0,X\in{}R^{N\times{}M}$$

With:

$$f_i\left(X\right)=x_{*,i}-F^{-1}\left[k\left(\omega{}\right)+\sum_{j=1}^Mq\left(\omega{}\right)\dot{}F\left[H^{\left(j\right)}\left(x_{*,j}\right)\right]\right]=0$$

Where $F$ is the Fourier transform and $H$ a nonlinear function, $k$ and $q$ are constants. $H$ has the following form:

$$H\left(x\right)=r{\left(\frac{x}{s}\right)}^{\alpha{}},\
s,\ r,\ x\in{}R^N$$

The highest $\alpha$ is, the harder it is to find a solution. I tried both simulated annealing and Jacobian-Free Newton Krylov numerical methods, without much success. I believe an issue here might be the way I state the problem.

Should JFKN or SA compute a perturbation/update on each matrix element separately? It does not sound reasonable to me, as I am dealing with timeseries, each column of $X$ being a timeseries of length $N$; I guess perturbation should be flexible yet consistent from one element to the next. Ha! A bit of insight or advice would be very much appreciated here!

Thank you