Take the 2-minute tour ×
Mathematics Stack Exchange is a question and answer site for people studying math at any level and professionals in related fields. It's 100% free, no registration required.

First off, long time reader, first time poster. Thanks in advanced for all the help this site has offered!

So the question!

I have two matrices in the form of the variables A and b that are inputs to my matlab function (posted below). I would like to calculate the significant figures used for the matrix inverse operation (matrix division) from the result A, b matrices. However, I have no idea where to start (matlab or mathmatically) to go about this approach. Help?

More context, using a square linear system (Ax=b) and I am seeing if its singular or nonsingular and trying to find a solution(s).

% x = answer
% y = 0 if no solution, 1 if nonsingular, 2 if many solutions
% z = p is number of sig figs
%
function [ x, y, z ] = squareLinSysSolv(A, b)


if det(A) == 0
    % Matrix is singular and therefor many solutions
    x = A\b;
    y = 0; % Used as place holder to compile
    z = 5; % Used as place holder to compile
elseif det(A) ~= 0
    % Matrix does not equal to zero (perhaps a number very close to it or
    % far from it) and therefor has a unique solution.
    x = A\b;
    y = 1; % Used as place holder to compile
    z = 5; % Used as place holder to compile
end
end

Edit: To make it a bit clear, z should be some integer which approximates (ceiling or floor value) a decimal number of significant figures at which A\b were calculated at.

Test Cases: Test/Spec sheet of what is expected. Both A and b are matrices and the result should be something like so.

A =
    1.5000    2.3000    7.9000
    6.1000    3.2000   13.0000
   13.0000   21.0000   76.0000

b =
     1
     3
     5
>> [x,y,z] = squareLinSysSolv(A,b)
% the result of x = A\b
x =

    0.8580
    3.0118
   -0.9132
% determinant is not equal to zero
y =

     1
% Amount of sig figs/precision in calculation
z =

     15
share|improve this question

1 Answer 1

There are a number of ways to do this. At its core, a floating point processor will compute an answer to a finite number of digits of precision. In most PCs, a standard double representation (which is what MATLAB would use), would give about 16 digits.

For instance, try this code and inspect the result:

clear all
eps
1-(1+eps/2)

The value eps is the machine epsilon -- essentially the smallest number that can be added to 1 and get an answer that is not 1. This value is around $2.2 \times 10^{-16}$, so that's where the 16 digits of precision comes from.

Now, when you compute A\b, MATLAB is calling one of many possible algorithms to compute the solution to $Ax=b$. Assuming that there is some error, we want to think about whether we can bound that error. A vector norm is a good way of looking at that:

$$b-x = e \implies |A^{-1}(b-x)| = |A^{-1}e|.$$

Computing the relative error, we may instead want to look at

$$\frac{||A^{-1}b||/||A^{-1}e||}{||b||/||e||}.$$

Using the properties of matrix and vector norms, we'll define a function of $A$, called the condition number:

$$\textrm{cond}(A) = ||A^{-1}||\cdot ||A||.$$

The matrix norm is very difficult to compute exactly, but it can be bounded:

$$ ||A|| \le \left(\sum_{i,j} a_{i,j}^2\right)^{1/2}. $$

What does this mean? Well, if you compute the matrix norm in this way, you can get a bound on the relative error of the solution, which should be enough for you to estimate the number of "significant digits," better known as "digits of precision." Indeed, the condition number is the upper bound of this relative error! So by computing the condition number, you know that your error will be no worse than that.

Finally, in MATLAB, to compute the condition number, type:

cond(A)
share|improve this answer
    
Thanks for the response and detailed explanation! I have a some-what good idea of what you are saying get I need the condition number, which will give me the a good enough estimate to the digits of precision (I like your wording better) However, when it comes to finding the conditional number, I am lost. I added a test case to my original post. When computing z I need to get a number close to 15 (depending on how liberal I am with rounding up to 16 or down to 14, etc). Doing the conditional math on the value of x keeps returning 1. I may be confused on the upper bound case. –  Adam Oct 25 '12 at 20:22
    
Think about this: if you have a true value $x$, and an approximately true value $x+\epsilon$, then the relative error $|x+\epsilon|/|x|$ tends to 1 as $\epsilon$ tends to zero. Therefore, the relative error is a good measure of how accurate the estimate is, and consequently, how many digits of precision difference the approximation will have -- e.g., if $x = 15$, look at the relative error estimates of $x_1 = 15.01$ and $x_2 = 15.0000001$. If $\textrm{cond}(A) \approx 1$, then your estimate x = A\b is basically as close to accurate as possible in double precision arithmetic. –  Arkamis Oct 25 '12 at 20:34

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.