GoHokies
  • Member for 6 years, 5 months
  • Last seen more than a week ago
"Cookbook" about iterative linear solvers and preconditioners
15 votes

Have a look at Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods (Barrett et al.). You can find it here. Here's why I'm recommending this over other references: the ...

View answer
Smallest eigenvalue without inverse
Accepted answer
13 votes

Compute the largest-magnitude eigenvalue $\lambda_{max}$ of $A$ (with, say, eigs('lm')). Then compute the largest magnitude (negative) eigenvalue $\hat{\lambda}_{max}$ of $M = A - \lambda_{max}I$ (...

View answer
Why is it not computationally possible to accurately predict the weather that would occur after 14 days?
12 votes

Further to Chris' answer: Yes, weather (or the equations describing it) is extremely sensitive to the initial conditions. The fact that the weather system contains phenomena at pretty much all time ...

View answer
Convexity of Sum of $k$-smallest Eigenvalue
Accepted answer
9 votes

Given $A \in {\bf S}^n$ (a positive definite matrix) with eigenvalues $\lambda_1 \leq \lambda_2 \leq \ldots \leq \lambda_n $, then: $\displaystyle f_k(A)=\sum_{i=1}^{k} \lambda_i$ is concave. Why? $$...

View answer
Eigenvectors of a small norm adjustment
7 votes

There exist special techniques for updating the eigen-decomposition of time-dependent covariance matrices. Given a "prior" eigenvalue decomposition (say at some initial time $t^0$), these recursive ...

View answer
How to start using LAPACK in c++?
7 votes

Here's another answer in the same vein as the above. You should look into the Armadillo C++ linear algebra library. Pros: The function syntax is high-level (similar to that of MATLAB). So no ...

View answer
Can ARPACK exploit hermiticity when diagonalising a complex matrix?
Accepted answer
6 votes

No, there is no specialized ARPACK routine for complex Hermitian matrices. The ARPACK authors recommend using the znaupd routine for both Hermitian and non-Hermitian problems: https://www.caam.rice....

View answer
Finite Elements: using preconditioned conjugate gradients with incomplete cholesky decomposition
Accepted answer
6 votes

While preconditioned CG with incomplete Cholesky (ICC) is reasonably straightforward to formulate mathematically, writing an efficient implementation is a non-trivial matter. Here's some of the ...

View answer
I've developed a derivative-free optimization method, looking for comments
6 votes

I would like to hear comments from users that have some practical models (e.g. black-box hyperparameter optimization) which are still needed to be solved acceptably - whether this method works or not ...

View answer
pointwise vs. continuous observations in PDE inverse problem
Accepted answer
6 votes

The measurements of this field are often spotty and missing chunks; why interpolate to get a continuous field of dubious fidelity if that can be avoided? You're perfectly right - most of the time, ...

View answer
How can i get gauss-lobatto points on a quadrilateral?
Accepted answer
6 votes

For domains that are logically square or cubic (like your quadrilateral), you can use the tensor product (dimension-by-dimension) approach. That is, generate the 2D Gauss-Lobatto point matrix as the ...

View answer
2d Euler manufactured solutions
Accepted answer
5 votes

This 2004 paper by Roy et. al in Int. J. Numer. Meth. Fluids 2004; 44:599–620 (DOI: 10.1002/d.660) should contain exactly what you're looking for: ftp://ftp.demec.ufpr.br/CFD/bibliografia/...

View answer
When training a neural network, why choose Adam over L-BGFS for the optimizer?
Accepted answer
4 votes

This topic has been discussed at some length on Cross Validated (aka stats.stackexchange) and Reddit: Why is Newton's method not widely used in machine learning? (see in particular Nick Alger's ...

View answer
Finite volume software packages
4 votes

Have a look at Pyclaw. This library has been around for quite a while and is fairly robust. It offers: Implementations of several Godunov-type methods and Riemann solvers in 1/2/3D. Adaptive mesh ...

View answer
adjoint method for reaction-diffusion problem
Accepted answer
4 votes

Some remarks: Notation: $[t_1,t_2]$: the time interval; $\Omega$: the spatial domain; $\bar{u}^i(x)$: the known tumour profile at $t_i$; $\left\| \cdot \right\|_\Omega$: a suitably chosen norm on $\...

View answer
Generate a set of orthogonal vectors to a given vector
4 votes

How about doing an SVD to get the null space of the $1 \times d$ matrix $A = [\bar{v}]$ with $\bar{v} := v^T/\|v\|$? The null function in MATLAB does exactly this: d = 4; v = rand(d,1); vbar = v.'/...

View answer
Is it necessary to invert precondition matrix for iterative solver?
Accepted answer
4 votes

OK - so your original equation is $$Ax = b$$ Say you've come up with a good preconditioner for $A$, call this $M$. Also, say you have pre-computed an LU-decomposition for this $M$, i.e. $$M = L_m ...

View answer
Finite difference recursion and higher order
Accepted answer
4 votes

Short answer: yes (in exact arithmetic). You'll have to use the centered difference formula evaluated at $x \pm \frac{1}{2}\delta x$, like this: $$ u_x = \frac{u(x + \frac{1}{2}\delta x) - u(x - \...

View answer
MATLAB : Does the qr algorithm and the DGEMM used in MATLAB take into account if the input matrix is tridigonal and optimize accordingly?
Accepted answer
4 votes

Here's a short Matlab benchmark. I've used a sparse tridiagonal matrix. function benchmark_qr maxit = 25; time = zeros(maxit,2); for iter = 1:maxit mdim = 500 * iter; A = gallery('...

View answer
Appropriate iterative linear solver for an eigenvalue problem
Accepted answer
4 votes

If your matrices are large, why not use a library like ARPACK? The shift-and-invert mode of ARPACK will help you calculate the eigenvalues close to $\sigma$. There are interfaces to ARPACK for most ...

View answer
Matrix vector multiplication performance
3 votes

An addition to the other two (excellent) answers: Understanding the source code of something like OpenBLAS can be a daunting task. As an alternative starting point, you could read through How to ...

View answer
A 95% minimal rectangle problem
Accepted answer
3 votes

Your problem reads Given $N$ points in, say, two dimensions, find an axes-parallel rectangle of minimal area that encloses at least 95% of the points. This is essentially equivalent to the ...

View answer
Preconditioning ARPACK eigenvalue solver
Accepted answer
3 votes

I've summarized the comments thread of your original question into an answer. Here are a few things that you can try: Increase the number of your Arnoldi vectors (NCV) generated at each iteration. ...

View answer
proper derivation of a functional for a time dependent parameter estimation problem
Accepted answer
2 votes

You want to solve $$ \arg \min_{D(x),k(x)} \mathcal{G}(u;D,k) := \frac{1}{2} \int_{\Omega} \left(u(t=T) - {\bar{u}}(T) \right)^2 d\Omega + \frac{1}{2}\int_{\Omega} D^2 d\Omega + \frac{1}{2}\int_{\...

View answer
Precision loss in Matrix-Vector product when applying Finite-Difference scheme
2 votes

If I understood you correctly, you are using the IDA solver from the LLNL SUNDIALS suite. If you have access to the source code that computes the residual of your DAE system, you may then consider ...

View answer
Implementing odespy for system of PDEs
Accepted answer
2 votes

Odespy needs from you a 1st order system of ODEs, namely something of the form \begin{align} \frac{d\mathbf{y}}{dt} = \mathbf{f}(\mathbf{y},t), \tag{1} \end{align} where, for your particular problem, ...

View answer
Numerical integration using RKF7(8) - different results
2 votes

One remark (seeing that the two answers given thus far come with quite a few comments, I've decided to post this as a separate answer instead): After you've solved the problem with your ODE ...

View answer
Sparse Matrix Reordering
2 votes

You can avoid the explicit computation of $A^{-1}$ in the Newton iteration for the sign of $A$. Refer to chapter 5 from this book by Higham, more specifically equation 5.22.

View answer
Automatic Differentiation - reverse accumulation of linear system solve
Accepted answer
1 votes

Starting with $b=Ax$ the forward differentiation gives $\dot b=\dot A x+A\dot x$, or, with $\dot y = \dot A x$: $$\dot x = -A^{-1} \dot y + A^{-1} \dot b$$ Let's write this in terms of the local ...

View answer
Why do equi-spaced points behave badly?
1 votes

It's good to be aware of Floater-Hormann interpolants when you have to (or want to) work with equidistant points $\{x_i\}_{i=1\ldots n}$. Given the integer $d$ with $0 \le d \le n$, let $p_i$ be the ...

View answer