Introduce the vector $y:=-A^{-1}Gx$ and solve the large coupled system $Ay+Gx=0$, $G^Ty=-b$ for $(y,x)$ simultaneously, using an iterative method. If $A$ is symmetric (as seems likely though you don't ...

No such standards exist, as reliable error estimates often cost much more than the approximate calculations. Basically there are four kinds of error estimates: (i) Theoretical analyses proving that ...

As $\det(kA)=k^n\det A$, the determinant can be made arbitrarily large or small by simple rescaling (which doesn't change the condition number). Especially in high dimensions, even scaling by an ...

Matlab interprets sequences of multiplications and/or divisions from left to right. Hence $A*B*C*v$ is much more expensive than $A*(B*(C*v))$, as you have two matrix products and one matrix-vecor ...

It depends a lot on the size of your matrix, in the large-scale case also on whether it is sparse, and on the accuracy you want to achieve. If your matrix is too large to allow a single ...

In low dimensions, a well implemented BFGS method is generally both faster and more robust than CG, especially if the function is not very far from a quadratic. Neither BFGS nor CG need any assumption ...

In Gerschgorin's theorem, the diagonal entries $A_{ii}$ of the matrix are the eigenvalue estimates, and the radii $r_i$ of the Gerschgorin disks are corresponding error bounds. Thus $\min_i A_{ii}-r_i$...

As your matrix is independent of $u$ the result is a matrix exponential times the intial vector. The standard discussion of relevent method can be found from http://scholar.google.at by searching for '...

The Cholesky factorization $C=R^TR$ leads to a Cholesky-like factorization of the inverse $C^{-1}=SS^T$ with the upper triangular matrix $S=R^{-1}$. In practice, is best to keep the inverse factored....

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 formula $$\mathrm{logsum}(x,y)=\max(x,y)+\mathrm{log1p}(\exp(-\operatorname{abs}(x-y))$$ should be numerically stable. It generalizes to a numerically stable computation of \log \sum_i e^{x_i} ...

Simple functions like Rosenbrock's are used to debug and pre-test newly written algorithms: They are fast to implement and to execute, and a method that cannot solve the standard problems well is ...

search = attempt to find a feasible point that satisfies all constraints (and for optimization a better point than found so far), generally using function values only. local search: improving a ...

CVXOPT only solves (smooth and nonsmooth) convex problems, giving access to several third party convex solvers with guaranteed state of the art worst case complexity. You may pose linear, convex ...

The standard way of doing it is to extract from the expression for $f(x)$ an exponential prefactor, transform that to $e^{-x^2}$, and then use Gaussian quadrature rules (or Gauss Kronrod) with this as ...

Adding squared penalty terms to get rid of constraints is a simple approach giving an accuracy of order 1/penalty factor only. Hence it is not recommended for high accuracy unless you let the penalty ...

Some references on rounding error analysis of Krylov methods: Gutknecht, Martin H., and Zdenvek Strakos. "Accuracy of two three-term and three two-term recurrences for Krylov space solvers." SIAM ...

To get a dense positive definite matrix with condition number $c$ cheaply, pick a diagonal matrix $D$ whose diagonal consists of numbers from $[1,c]$ (which will be the eigenvalues), with $1$ and $c$ ...

Checking feasibility of an LP and solving an LP are essentially equivalent problems, as one can be transformed into the other by standard methods changing the complexity by a constant factor only. ...

Their observation ''the form of the algorithm may needlessly introduce some numerical imprecision'' is correct. But their explanation ''This arises because, in eqn (3.14), a small term ($O(δt^2)$) is ...

Usually one doubles the initial step until the Goldstein condition is violated or (in a feasible point method) the boundary is reached. Then one has a bracket. (If no such step exists, the objective ...

TQL cannot be parallelized. The standard parallel algorithm is that of Cuppen: JJM Cuppen, A divide and conquer method for the symmetric tridiagonal eigenproblem, 1980. http://www.springerlink.com/...

If $m\gg n$, as your question indicates, you can save some work by first picking an index set $I$ of $p\approx 5n$ (say) random rows and using the orthogonal factorization $A_{I:}^T=QR$. (The QR-...

Strictly speaking, the problem of computing multiplicities is ill-posed, as arbitrarily small perturbations may change the multiplicities (usually reducing them to 1). However, to some approximation, ...

A sparse direct solver knows the matrix, and hence its dimensions and its sparsity pattern. Of course it doesn't know the dimension of the problem dimension before discretization. However, the ...

This is far from generalized SVD. If B is a positive matrix, you can use my package BIRSVD http://www.mat.univie.ac.at/~neum/software/birsvd/ The paper http://www.mat.univie.ac.at/~neum/software/...

See L.M. Rios and N.V. Sahinidis, Derivative-free optimization: A review of algorithms and comparison of software implementations for a very useful recent comparison of solvers. DOI: 10.1007/...

Yes. Compute the $QR$ factorization and take $L=R^T$; rescale the rows of $R$ if necessary (by changing some of their signs) to make the sign of the diagonal nonnegative (as the Cholesky factor is ...

In theory, yes. In practice, rounding errors will usually result in (initially slow) convergence to $u_1$. At essentially the same cost one can run the Lanczos algorithm, which will have much faster ...