26

Yes, the condition number always matters in floating-point arithmetic, whether you choose to solve your system with an iterative or direct method. The relative accuracy of an approximate solution to $Ax = b$ obtained from LU factorization with pivoting is $O(\kappa(A) \cdot \varepsilon)$, where $\varepsilon$ is the smallest floating point number such that $1 ...


13

You can use additive $$ P_a^{-1} x = (B^T B)^{-1} x + (C^T C)^{-1} x, $$ multiplicative $$ P_m^{-1} x = (B^T B)^{-1} x + (C^T C)^{-1} \Big(x - A (B^T B)^{-1} x \Big), $$ or symmetric multiplicative. Methods of this class are available in PETSc using PCCOMPOSITE in PETSc. For example, petsc/src/ksp/ksp/examples/tutorials$ ./ex2 -m 100 -n 100 -...


13

When you use ZGELSS to sovle this problem, you're using the truncated singular value decomposition to regularize this extremely ill-conditioned problem. it's important to understand that this library routine is not attempting to find a least squares solution to $Ax=b$, but rather it is attempting to balance finding a solution that minimizes $\| x \|$ ...


12

Though it is a relatively rare situation when you actually have to calculate an inverse of the matrix, not all techniques were created equally. I would use the term badly-conditioned instead of ill-conditioned. For badly conditioned matrices, you might opt in the SVD-route to calculate the inverse: $$ A=U\Sigma V^H \implies A^{-1}=V\Sigma^{-1}U^H. $$ If ...


11

MATLAB has a couple of "exact" functions for this, cond and rcond, with the latter returning a reciprocal of the condition number. Matlab approximate function condest is more fully described below. Often estimates of the condition number are generated as by-products of the solution of a linear system for the matrix, so you might be able to piggyback the ...


9

Solving a (linear) PDE consists in discretizing the equation to yield a linear system, which is then solved by a linear solver whose convergence (rate) depends on the condition number of the matrix. Scaling the variables often reduces this condition number, thus improving convergence. (This basically amounts to applying a diagonal preconditioner, see ...


9

PDEs of which the solutions have sharp boundaries pose problems that go beyond being able to represent the solution in floating point. This is especially true when solutions have a certain physical meaning, e.g., a density (that per se cannnot be smaller than 0). Consider, for example, $$ -\varepsilon \Delta u + u = 0 \text{ on } \Omega,\\ u = 1 \text{ on } \...


9

Jed Brown has already pointed this out in the comments to the question, but there is really not very much you can do in usual double precision if your condition number is large: in most cases, you will likely not get a single digit of accuracy in your solution and, worse, you can't even tell because you can't accurately evaluate the residual corresponding to ...


9

The condition number of sum $s(x) = \sum_{j=1}^n x_j$ is given by $$ \kappa(x) = \frac{\sum_{j=1}^n |x_j|}{|\sum_{j=1}^n x_j|} = \frac{s(|x|)}{|s(x)|}$$ and reflects the sums sensitivity to small changes in the input. Specifically, we have $$ \underset{\epsilon \rightarrow 0_+}{\lim}\sup \left\{ \frac{1}{\epsilon} \left|\frac{s(x+\Delta x) - s(x)}{s(x)} \...


8

Added after my initial answer: It appears to me now that the author of the referenced paper is giving condition numbers (apparently 2-norm condition numbers but possibly infinity-norm condition numbers) in the table while giving maximum absolute errors rather than norm relative errors or maximum elementwise relative errors (these are all different measures.) ...


8

No. A reordering of a matrix $\boldsymbol{A}$ is equivalent to conjugation by a permutation matrix $\boldsymbol{P}$. In other words, the reordered matrix can be written as $\boldsymbol{A}_r = \boldsymbol{P}\boldsymbol{A}\boldsymbol{P}^*$ and $\boldsymbol{A}_r^{-1} = \boldsymbol{P}\boldsymbol{A}^{-1}\boldsymbol{P}^*$. Since the spectral norm is unitarily ...


7

The simplest/fastest way to solve ill-conditioned problems is to increase precision of computations (by brute force). Another (yet not always possible) way is to re-formulate your problem. You might need to use quadruple precision (34 decimal digits). Even though 20 digits will be lost in a course (because of condition number) you will still get 14 correct ...


6

No, Jacobi only ever corrects relative scales. It does nothing for "smooth" ill-conditioning, such as the $\kappa(A) \in O(h^{-2})$ asymptotics for second order elliptic problems. If you are using a Krylov method, the global scale is automatically corrected, but with a stationary iteration, the (constant) scaling is needed somehow (could just be in the ...


6

Let's show that we cannot bound the condition number of $X^T A X$ by using only the condition number of $A$ and the norm of $X$. Let $A=I$, so its condition number is exactly $1$. Let $X$ consist of an invertible diagonal block with $p$ rows of zeros padded at the bottom: $$ X = \begin{pmatrix} D \\ 0 \end{pmatrix} $$ Now $X^T A X = D^2$, and its ...


6

The condition number measures the relative distance from singularity of a matrix $A$: that is, $$ \min_{\text{$X$ is singular}} \frac{\|A-X\|}{\|A\|} = \frac{1}{\kappa(A)} $$ (the norm here is the Euclidean / induced / spectral norm --- i.e., $\|A\|=\sigma_1(A)$). This property follows from the Eckart-Young theorem. Your example with a small $\|A\|$ ...


6

The condition number for the stiffness matrix of any method (finite elements, finite volumes, finite differences) applied to a second order differential operator always grows as ${\cal O}(h^{-2})$ where the exponent equals the order of the differential operator. As a consequence, if you just make the mesh fine enough, you can make the condition number ...


6

The inequality $\frac{\alpha + \lambda_{max}(A)}{\alpha+\lambda_{min}(A)} > \frac{\lambda_{max}(A)}{\lambda_{min}(A)}$ doesn't hold because $\lambda_{max}(A) > 1 > \lambda_{min}(A) > 0$ and $\alpha > 0$. In fact, the reverse inequality holds in this case. Hence, there is no contradiction here. EDIT: From the comments, I realized that the ...


5

No, a different unit will not alter the condition of the system. Say your FEM system is $$ Au + \beta Mu = 0. \quad (*) $$ Then the parameter $\beta$, here something like $b/a$ from your example, will depend on the units in a way that only allows you to add "$A$" and "$M$" in terms of units. A rescaling will then mean a multiplication of $(*)$ by a constant ...


5

Your question really doesn't admit a simple answer- we need to know more specifics about your problem to provide a useful answer. In general, iterative methods can be faster than direct factorization for large sparse systems of equations if the system is reasonably well conditioned or if it is badly conditioned but you have a good preconditioner or if you ...


5

The accuracy of Gaussian elimination is not bounded in terms of condition number. There is an example of a well-conditioned matrix in Trefethen and Bau (and probably elsewhere) for which Gaussian Elimination is exponentially unstable with respect to matrix size. Now for iterative methods: CG finds the best approximation over the subspace in the energy norm ...


5

Of course you can. The condition number gives you an upper bound for the error. Suppose you want to solve $Ax=b$ where $A$ is a matrix with large condition number. The error depends on the structure of $A$ and on $b$. As an example, if $A$ is a $2 \times 2$ diagonal matrix, with elements $1$ and $10^{-20}$, solving $Ax=b$ will still be accurate, despite ...


5

This algorithm is most useful for two situations, which are related to each other in practice: You don't know the matrix entries explicitly, but instead can only compute matrix-vector products with the matrix (often called "matrix-free") You want the 1-norm for the inverse of a matrix. The inverse of a sparse matrix is typically dense, and since sparse ...


5

Ill conditioning is a property of the system of equations rather than the algorithm used to solve the system of equations. Using a bad algorithm can certainly make the situation worse, but you're already in trouble when you try to solve an ill-conditioned system of equations with A coefficients or right-hand side b with even tiny errors even if you use ...


5

There is no simple fix. For an ill-conditioned matrix $A$, the harm (loss of precision) is already done the moment you wrote those numbers in a numpy array, because that tiny $10^{-16}$ perturbation from the exact non-representable values is already harmful. You could increase your working precision; but at that point the question is if your matrix entries $...


4

The Laplacian is SPD so the singular values are the same as the eigenvalues (and the eigenvectors are orthogonal). This article explains the eigenvalues and eigenvectors of the discrete Laplacian obtained using centered differences. http://en.wikipedia.org/wiki/Kronecker_sum_of_discrete_Laplacians


3

Please notice that you assumptions do not preclude singular matrices. A specific example is \begin{equation} C = \begin{bmatrix} I & I \\ I & I \end{bmatrix} \end{equation} where $I$ denotes the identity matrix of dimension $I$. This example also shows that the condition that $B$ is positive definite is not necessarily useful. I appreciate the fact ...


3

You can get a rough estimation of the condition number using the Gershgorin circle theorem. This article in Wikipedia has a nice explanation: http://en.wikipedia.org/wiki/Gershgorin_circle_theorem For a complex matrix $A$ of size $n\times n$, its n Gershgorin circles are drawn in the complex plane, the center of the circles are $a_{ii}$, the radius of the ...


3

Estimates can be obtained from Krylov subspace methods, such as GMRES, and PETSc has functionality for that. For a symmetric matrix, the 2-norm condition number is the spectral condition number, and there exist algorithms to estimate the 2-norm condition number, so if you really do want the condition number for your unpreconditioned matrix, you can use ...


3

(This just addresses the first part of your question). The subtlety here is in your definition of the condition number -- as defined here, it is a ratio of relative error in the output to relative error in the input (ie these notes). As you've discovered, near $x=a$, the relative error of the input is order $\epsilon$, while the relative error of the ...


3

Assuming your jacobian isn't too large, it's a good idea to solve for your newton direction by QR factorization when your jacobian has a large condition number. It is numerically stable since it breaks the problem into two separate systems with the minimum error amplification factor allowed by the original problem. See here for details. This enables you to ...


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