17

Forming the Schur complement Suppose that you have permuted and partitioned your matrix into the form $$A=\left(\begin{array}{cc}A_{11} & A_{12} \\ A_{21} & A_{22}\end{array}\right),$$ such that $A_{22}$ contains your degrees of freedom of interest and is much smaller than $A_{11}$, then one can form the Schur complement $$ S_{22} := A_{22} - A_{...


15

Although mixed-integer linear programming (MILP) is indeed NP-complete, there are solvable (nontrivial) instances of mixed-integer linear programming. NP-complete means that mixed integer linear programming is: a) solvable in polynomial time with a nondeterministic Turing machine (the NP part) b) polynomial time reducible to 3-SAT (the complete part; for ...


15

If your graph is undirected (as I suspect), the matrix is symmetric, and you cannot do anything better than the Lanczsos algorithm (with selective reorthogonalization if necessary for stability). As the full spectrum consists of 100000 numbers, I giess you are mainly interested in the spectral density. To get an approximate spectral density, take the ...


13

A fairly simple method would be to choose a basis in function space and convert the integral transformation to a matrix. Then you can just invert the matrix. Mathematically, here's how that works: you need some set of orthonormal basis functions $T_i(x)$. (You can get away without them being normalized too, but it's easier to explain this way.) Orthonormal ...


13

As others have pointed out, this is difficult to do with a direct solver. That said, it isn't that hard to do with iterative solvers. To this end, note that most iterative solvers in one way or another minimize the error with respect to some norm. Oftentimes, this norm is either induced by the matrix itself, but sometimes it is also just the l2 vector norm. ...


12

For the sake of notation, let's suppose that $f: \mathbb{R}^{n} \rightarrow \mathbb{R}^{n}$ (i.e., it's a vector-valued function that takes a vector as input and outputs a vector of the same size). There are two concerns: computational cost and numerical accuracy. Calculating the derivative $\mathrm{D}f(x)$ (the Jacobian matrix, $J(x)$, or $(\nabla f(x))^{T}...


9

There is a prime missing both in your code and in the expression in your question. In the original paper, the expression is: $$ \log(x/x_0) + (x-x_0) \frac{\sum^\prime_{0\leq j\leq n}p_j T_j^*(x/6)}{\sum^\prime_{0\leq j\leq n}q_j T_j^*(x/6)}. $$ This primed sum is very common and standard when dealing with Chebyshev series: it means that the first term of ...


8

The usual approach to testing optimization algorithms is to compare how many function evaluations they need to find the minimum (or get within a fixed tolerance $\varepsilon$ of the minimum). This is easily implemented in any code, and it does not matter in that case how fast your computer is, or what the resolution of the tic/toc mechanism is. There are a ...


6

The model reduction approach Since Paul asked, I'll talk about what happens if you use projection-based model reduction methods on this problem. Suppose that you could come up with a projector $\mathbf{P}$ such that the range of $\mathbf{P}$, denoted $\mathcal{R}(\mathbf{P})$, contains the solution to your linear system $\mathbf{Ax} = \mathbf{b}$, and has ...


6

One way to get the $d_k$ is to expand each Chebyshev polynomial in powers, and then take linear combinations. However, using the Chebyshev representation itself (as detailed below) is numerically far more stable than using the power series expansion. You can evaluate a Chebyshev sum at any pareticular $x$ in a Horner-like fashion by using the 3-term ...


6

Arnold Neumaier's answer is discussed in more detail in section 3.2 of the paper "Approximating Spectral Densities of Large Matrices" by Lin Lin, Yousef Saad and Chao Yang (2016). Some other methods are also discussed but the numerical analysis at the end of the paper shows that the Lanczos method outperforms these alternatives.


5

It may look like the error is getting larger, but I think you are confusing multiple issues here. The first is that you are looking at finite approximating Taylor polynomials $p_3(x), p_5(x), p_7(x)$ evaluated at $x=\pi/2$. Taylor's theorem says that for analytic functions (which the sine is), $p_N(\pi/2)\rightarrow 1$ as $N\to \infty$. Your results appear ...


5

The "fast" version is (damped/over-relaxed) Jacobi instead of SOR, which is a multiplicative method. Even with damping factor of 1.0, Jacobi is not guaranteed to converge, as you would see if you applied it to a 1D problem, or added extreme anisotropy to your multi-dimensional problem. Your SOR ("slow") implementation is slow because it is written in pure ...


5

You could try using a library that implements the Fast Multipole Method (FMM), which should drastically reduce the amount of memory you need and will decrease the complexity of matrix-vector products from $\mathcal{O}(N^{2})$ to $\mathcal{O}(N)$. It is difficult to implement, but there should be some libraries out there. Another algorithm for N-body ...


5

In order to inconvenience as many people as possible, long ago, mathematicians and physicists decided to use two different conventions on whether $\theta$ or $\phi$ is the latitude angle. Greengard's notes use the physicists' convention that $\theta$ is latitude and $\phi$ is longitude, whereas Scipy uses $\theta$ for longitude and $\phi$ for latitude, so ...


5

You can write a curve as a parametric equation ($x(t)$, $y(t)$) for a range say $0 < t < 1$. Then you can have a Fourier decomposition of both $x(t)$ and $y(t)$. A rectangle can be written as such a parametric equation. However, you want to decompose a given function as a sum of functions from a specific class. There is a whole field devoted to this. ...


5

I hope I understood the question correctly. They try to compute exactly the same thing, so they really are equivalent. I'll use Chebyshev polynomials because they are easy to analyze. Given a function $f(x)$ on $[-1,1]$, the spectral interpolant is the truncation of $$ \begin{aligned} f(x) &= \sum_{n\geq0} \bar a_n T_n(x), \\ \bar a_n &= \frac{1+[n&...


5

I believe that for odd polynomial degrees you can get superconvergence by one order at the jump. For even degrees I don't think that's the case. For the odd case, we consider a smooth function $u:\mathbb{R}\to\mathbb{R}$. Consider two adjacent intervals $I^-=[x_{i-1},x_i]$ and $I^+=[x_i,x_{i+1}]$ both of length $h$, and denote the $L^2$ projection of $u$ on ...


5

I think this is a really cool question, which might have a really cool answer, if someone was willing to think about it hard enough to publish on. But when we're talking smooth functions which need to be evaluated quickly, I reach for B-splines. Create the pairs $\{x_j, f(x_j)\}_{j=1}^{n}$, and use (say) a cubic B-spline (here's an implementation for the ...


4

If you're ok with thinking about things that are not eigenvalues but functions that in some sense still tell you something about the spectrum, then I think you should check out some of the work by Mark Embree at Rice University.


4

The long answer is...sort of. You can re-arrange your system of equations such that the farthest right $k$ columns are the variables which you wish to solve for. Step 1: Perform Gaussian Elimination so that the matrix is upper triangular. Step 2: solve by back substitution for only the first (last) $k$ variables which you are interested in This will save ...


4

The "right" answer strongly depends on what you need your approximant for. Do you really need the best approximation for some error bound? Or just a good approximation? Or just a good approximation in the minmax sense? Nick Trefethen recently gave a nice example where Remez approximation is a bad idea since it minimizes the maximum error irrespective of the ...


4

There are many different types of approximations (or "surrogate models") you could try. Some that come to mind are Kriging, MARS, and Radial Basis Functions. These types of surrogate models (as opposed to polynomial regression) can accommodate a wide range of functional relationships, but you might need to experiment a bit to find which works best for your ...


4

For simplicity's sake, let's just say that $a=0,b=1$. Then you want to estimate $$ \int_0^1 f(x) g(x) dx - \left(\int_0^1 f(x) dx \right) \left(\int_0^1 g(x) dx\right) $$ from above and below. Furthermore, let me consider two cases: If $g(x)$ has mean value zero, i.e., $\int_0^1 g(x) dx = 0$, then you know a priori that your approximation isn't a good one....


4

Let's use the example (approximate the square root function on $[0.25,1.0]$ with a quartic polynomial) to step through your calculations. I suspect that the code is going to work with only modest changes. Having chosen an initial set of six interpolation points: $$ x_k = 0.25, 0.4, 0.55, 0.7, 0.85, 1.0\;\; (k=0,\ldots,5) $$ we proceed to interpolate the ...


4

Randomized algorithms can accurately approximate your matrix if it has rapidly decaying singular values (i.e. a Gaussian sampling of the range is likely to pick up the action of the matrix). Since your matrix is symmetric, you could approximate a truncated eigenvalue decomposition using random linear algebra and simply take the square root of the eigenvalues....


4

Taylor series is not a good way to do this; it takes a lot of terms to get reasonable accuracy. I answered a similar question on stackoverflow a while ago. Here's the body of the answer: Here are some good slides on how to do power series approximations (NOT Taylor series though) of trig functions: http://www.research.scea.com/gdc2003/fast-math-...


3

I'd follow your suggestion of fitting a line $ax+b$ to your data. You don't even have to smooth it first -- fitting a line already takes care of it. There are, however, two questions: first, how do you fit the line? This is a question of the statistics of your noise. If noise is Gaussian then you'd use least squares. If your noise has large outliers then ...


3

Showing that one problem can be reduced to another to demonstrate that it is NP-hard is in general not a trivial matter. It means that you have to show that you can reformulate the problem as another one for which it is known that it is NP-hard and then obtain the solution of the original problem from the solution of the reformulated problem. The more ...


3

Write your function as a linear combination of harmonics and treat the coefficients as variables. This gives a semi-infinite linear program. The Matlab optimization toolbox has a routine linprog to solve linear programs (LPs) only; so you need to discretize the curvature constraint. Replace the curvature constraint by $N$ constraints evaluated at $N$ ...


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