22

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 is a matrix ill conditioned? It depends on the accuracy of the solution you are looking for, as much as "beauty is in the eye of the beholder"... May be your question should better rephrased as are there cheap and robust condition number estimators based on the $LU$ factorization? Assuming you are interested in the real general (dense, non symmetric) ...


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

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 computing a few such solutions will tell you what is going on.


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

Computing the condition number (even approximating it within a factor of 2) seems to have the same complexity as computing a factorization, though there are no theorems in this direction. From a sparse Cholesky factor $R$ of a symmetric positive definite matrix, or from a sparse $QR$ factorization (with implicit $Q$) of a general square matrix, one can ...


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


7

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


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

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

Preliminary observations If you have a basis $\{\phi_i(x)\}$ of size $n$, you must have at least $n$ quadrature points. You can see this by writing the inner product as $(u,v) = u^T B^T W B v$ where $W$ is a diagonal weight matrix and $B_{ij} = \phi_j(x_i)$ is the basis evaluation matrix. If there are fewer than $n$ quadrature points, $B$ has fewer than $n$ ...


5

It is nearly impossible to to tell if the your system is ill conditioned from just one result. Unless you have some foresight into the behavior of your system (i.e. know what the solution SHOULD be), there's not much you can say from a single solution. Having said this, you can gain more information if you solve more than one system with the same $A$. ...


5

Multiply the diagonal elements by a factor $q>1$ but close to 1. It will usually do the job. (I'd first try $q=1.0001$, but one can experiment with the number of zeros in this expression; e.g., use $q=1.01$ for very noisy data.) This is justified under certain conditions, as it is a specific form of regularization. For more on regularization, see my ...


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

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


4

The question to ask is why you need to invert the matrix. If a matrix is near-singular, it's true that you can define something like a pseudo-inverse in some stable way but it's nevertheless true that for near-singular matrices, solving linear systems in any way is unlikely to result in anything useful because the result is so strongly dependent on small ...


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