Search Results

Answering to the edit: I still have no idea how $S$ is defined in terms of $G$. But that's fine, it can just be a $D^2\times D^2$ matrix for now. Not all matrices can be written as $S = U \otimes V$. …

A Kronecker product is essentially a rank-1 matrix, once you change indices: if $A_{ijkl} = U_{ik}V_{jl}$ (where all indices vary from $0$ to $n-1$, to keep the notation simpler), then the matrix defi …

Since you speak about mathematical proofs, probably you don't want an arbitrary-precision library, but a verified computation library based on interval arithmetic. Maybe Arb? Or boost-interval? And ma …

[EDIT: not working in fact because the update of $Q$ cannot be merged efficiently, see comments] A simple eigendecomposition $A=QDQ^*$ should work, since it can be updated in $O(n^2)$.

It looks like you want $m - QQ^Tm$, is that correct?

No, the first line is not equivalent to what you claim in the text, I believe: that's one linear system containing a subtraction and not two separate linear systems. In detail, the matrix of this line …

Based on what I know of similar equations, I would expect that you can't beat this "successive substitution" (essentially a non-linear Gauss-Seidel) in case where it converges fast and only few iterat …

A variant of the argument in Carlo Beenakker's answer: if the $x_i$ are equispaced points with distance $h$ one from the next, then $$\frac{f(x_{i+1})-f(x_i)}{h} - f'(x_i) = O(h),$$ $$ \frac{f(x_{i+1} …

You are essentially using normal equations to solve the least-squares problem $\min \|W^{1/2}(Ax-b)\|_2$ resulting from IRLS. Normal equations are known not to be a backward stable algorithm. Use othe …

I suspect that using the expm1 function would give you a better result. Computing $e^x -1$ with the trivial formula in machine precision gives you only limited accuracy for small inputs: the fundament …

Turning my previous comments into an answer. What you have is an algebraic Riccati equation: indeed, setting $F=-\frac12 I$, you get $F^TX+XF + B = XAX$. Since $A$ and $B$ are positive definite, the …

With the traditional algorithms and complexity measures used in numerical linear algebra (dense real matrices, floating point computations, flop count as a complexity measure), they are both more or l …

Disclaimer 1: Treating these topics properly would require a quick course in numerical analysis. Disclaimer 2: If you are using any sane computer system, it's already going to have a library function …

At least in my field (numerical linear algebra), the current standard is that including the full source code is not mandatory for a publication. That said, there are many reasons why sharing your code …

Expand $(M3 - \frac13 \operatorname{Tr}(M3))^3 = 0$. This gives you simpler equations than the ones you are using. Note that this equation must hold in the triple point, because $M3 - \frac13 \operat …