Arnold Neumaier
  • Member for 9 years, 10 months
  • Last seen more than a week ago
Solving $(G^TA^{-1}G)x = b$ without inverting $A$
Accepted answer
22 votes

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 ...

View answer
Scientific standards for numerical errors
19 votes

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 ...

View answer
Does a tiny determinant imply ill-conditioning of a matrix?
19 votes

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 ...

View answer
MATLAB matrix multiplication (the best computational approach)
17 votes

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 ...

View answer
What is the fastest way to calculate the largest eigenvalue of a general matrix?
17 votes

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 ...

View answer
BFGS vs. Conjugate Gradient Method
16 votes

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 ...

View answer
Gershgorin Circle Theorem to estimate the eigenvalues
Accepted answer
15 votes

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$...

View answer
Algorithms for linear system of ODEs
Accepted answer
15 votes

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 '...

View answer
Dealing with the inverse of a positive definite symmetric (covariance) matrix?
Accepted answer
15 votes

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....

View answer
Approximate spectrum of a large matrix
15 votes

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 ...

View answer
Catastrophic cancellation in logsum
Accepted answer
13 votes

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} ...

View answer
Testing numerical optimization methods: Rosenbrock vs. real test functions
13 votes

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 ...

View answer
Meaning of search methods and optimization methods
Accepted answer
12 votes

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 ...

View answer
CVXOPT VS. OpenOpt
Accepted answer
11 votes

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 ...

View answer
How can I approximate an improper integral?
11 votes

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 ...

View answer
Nonlinear least squares with box constraints
11 votes

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 ...

View answer
Krylov subspace iterative methods in floating point arithmetic
Accepted answer
10 votes

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 ...

View answer
Generating Symmetric Positive Definite Matrices using indices
Accepted answer
10 votes

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$ ...

View answer
LP feasibility checking
9 votes

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. ...

View answer
In floating point arithmetic, why does numerical imprecision result from adding a small term to a difference of large terms?
Accepted answer
9 votes

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 ...

View answer
Initially Bracketing Minimum for Line Search
Accepted answer
9 votes

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 ...

View answer
Parallel algorithm for eigensystem of a tridiagonal matrix
9 votes

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/...

View answer
Null-space of a rectangular dense matrix
9 votes

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-...

View answer
How to detect the multiplicity for the eigenvalues?
Accepted answer
9 votes

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, ...

View answer
How does a Sparse Direct Solver know about dimensionality of a problem being solved?
Accepted answer
9 votes

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 ...

View answer
weighted SVD problem?
Accepted answer
9 votes

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/...

View answer
Finding a global minimum of a smooth, bounded, non-convex 2D function that is costly to evaluate
8 votes

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/...

View answer
Computation of Cholesky factor
8 votes

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 ...

View answer
What are the symptoms of ill-conditioning when using direct methods?
8 votes

The solution of an ill-conditioned system of equations with a matrix of norm 1 a random right hand side of norm 1 will have with high probability a norm of the order of the condition number. Thus ...

View answer
Eigenvectors with the Power Iteration
8 votes

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 ...

View answer
1
2 3 4 5
7