inv(A)*B*inv(A) or A\B/A, which is more accurate for a full rank A and half rank B

My matrices are

A = [0	0	0	0	0	0	0	0	0	1
0	0	0	0	0	0	0	0	1	0
0	0	0	0	0	0	0	2	0	0
0	0	0	0	0	0	6	0	0	0
0	0	0	0	0	24	0	0	0	0
1953050.78182336	390611.805629631	78122.6909803108	15624.6041672176	3124.93402773033	624.989444414569	124.998416658843	24.9997888874000	4.99997888869543	1
3515506.25066668	624981.527842487	109372.229170523	18749.6041663820	3124.94722207284	499.993666635371	74.9993666622000	9.99995777739086	1	0
5624833.75058238	874977.833364188	131247.229164674	18749.6833324371	2499.96833317685	299.997466648800	29.9998733321726	2	0	0
7874800.50027769	1049977.83331739	131247.783327059	14999.8099990611	1499.98733324400	119.999493328690	6	0	0	0
9449800.49985652	1049982.26661648	104998.669993428	8999.92399946400	599.997466643452	24	0	0	0	0];

B=[0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0
3515506.25066668	624981.527842487	109372.229170523	18749.6041663820	3124.94722207284	499.993666635371	74.9993666622000	9.99995777739086	1	0
5624833.75058238	874977.833364188	131247.229164674	18749.6833324371	2499.96833317685	299.997466648800	29.9998733321726	2	0	0
7874800.50027769	1049977.83331739	131247.783327059	14999.8099990611	1499.98733324400	119.999493328690	6	0	0	0
9449800.49985652	1049982.26661648	104998.669993428	8999.92399946400	599.997466643452	24	0	0	0	0
9449840.39954828	839989.359947423	62999.4679962480	3599.98479986071	120	0	0	0	0	0];

I tried both inv(A)*B*inv(A) or A\B/A, resulting in very different results for a gradient descent algorithm.

So, which is more accurate?

[Torsten]

Did you think about why you get different results ? Is your matrix A near to singular and badly conditioned ? What do you get with

cond(A)

?

Can you include the matrix here for testing ?

[Qiang Li]

Hi Torsten, here's cond(A)

K>> cond(A_2stepblock)

ans =

     1.439014477720320e+07

I updated the question with the matrices.

Answer 1

[John D'Errico]

Why would you want to compute the inverse of A TWICE ANYWAY? Even using both slash and backslash is arguably a bad idea, since again, you are doing way more work than you need.

Instead, decompose A. For example, you might decide to compute a QR factorization for A. Or perhaps an LU. Even then, the condition number of A is large enough that when you do two solves, it still inflates any noise in the system by the SQUARE of the condition number. And since you have given us only numbers to 14 significant digits, that means anything you get will probably be complete crap, no matter what you do.

A = [0	0	0	0	0	0	0	0	0	1
    0	0	0	0	0	0	0	0	1	0
    0	0	0	0	0	0	0	2	0	0
    0	0	0	0	0	0	6	0	0	0
    0	0	0	0	0	24	0	0	0	0
    1953050.78182336	390611.805629631	78122.6909803108	15624.6041672176	3124.93402773033	624.989444414569	124.998416658843	24.9997888874000	4.99997888869543	1
    3515506.25066668	624981.527842487	109372.229170523	18749.6041663820	3124.94722207284	499.993666635371	74.9993666622000	9.99995777739086	1	0
    5624833.75058238	874977.833364188	131247.229164674	18749.6833324371	2499.96833317685	299.997466648800	29.9998733321726	2	0	0
    7874800.50027769	1049977.83331739	131247.783327059	14999.8099990611	1499.98733324400	119.999493328690	6	0	0	0
    9449800.49985652	1049982.26661648	104998.669993428	8999.92399946400	599.997466643452	24	0	0	0	0];
cond(A)

ans = 1.4390e+07

cond(A)^2

ans = 2.0708e+14

Do you see that any noise in the least significant digits of B will potentially be amplified by the SQUARE of the condition number here, and that means, since you have only 14 significant digits in both A and B, the result may be meaningless?

B=[0	0	0	0	0	0	0	0	0	0
    0	0	0	0	0	0	0	0	0	0
    0	0	0	0	0	0	0	0	0	0
    0	0	0	0	0	0	0	0	0	0
    0	0	0	0	0	0	0	0	0	0
    3515506.25066668	624981.527842487	109372.229170523	18749.6041663820	3124.94722207284	499.993666635371	74.9993666622000	9.99995777739086	1	0
    5624833.75058238	874977.833364188	131247.229164674	18749.6833324371	2499.96833317685	299.997466648800	29.9998733321726	2	0	0
    7874800.50027769	1049977.83331739	131247.783327059	14999.8099990611	1499.98733324400	119.999493328690	6	0	0	0
    9449800.49985652	1049982.26661648	104998.669993428	8999.92399946400	599.997466643452	24	0	0	0	0
    9449840.39954828	839989.359947423	62999.4679962480	3599.98479986071	120	0	0	0	0	0];

Even so, we might try a QR or an LU. Simplest would be to do this:

Afac = decomposition(A)

Afac = 
  decomposition with properties:

    MatrixSize: [10 10]
          Type: 'lu'

  Show all properties

As you can see, it chose an LU by default.

Afac\B/Afac

ans = 10x10

-0.0001 -0.0001 -0.0001 -0.0001 -0.0000 0.0001 -0.0001 0.0001 -0.0001 0.0000 0.0013 0.0029 0.0027 0.0013 0.0003 -0.0013 0.0028 -0.0025 0.0011 -0.0002 -0.0097 -0.0215 -0.0202 -0.0096 -0.0020 0.0097 -0.0200 0.0170 -0.0070 0.0012 0.0323 0.0717 0.0672 0.0320 0.0067 -0.0323 0.0627 -0.0493 0.0184 -0.0027 -0.0403 -0.0896 -0.0840 -0.0400 -0.0083 0.0403 -0.0717 0.0504 -0.0160 0.0017 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

That requires MATLAB to perform only one operation in the form of a decomposition. As such, it will be considerably more efficient to perform. It will be a far better solution that using the inverse matrix almost always.

My real question is, can you improve the condition of A? If I look at the last 5 rows of A, it looks like a Vandermonde matrix. That would suggest re-scaling the polynomial variable to be on the order of 1, instead of going out as far as 5. All it is is a variable re-scaling, but it would improve the problem dramatically. For example...

Ahat = [A(1:5,:);A(6:10,:)*diag(5.^(-9:1:0))]

Ahat = 10x10

0 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 0 1.0000 0 0 0 0 0 0 0 0 2.0000 0 0 0 0 0 0 0 0 6.0000 0 0 0 0 0 0 0 0 24.0000 0 0 0 0 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.7999 1.6000 1.4000 1.2000 1.0000 0.8000 0.6000 0.4000 0.2000 0 2.8799 2.2399 1.6800 1.2000 0.8000 0.4800 0.2400 0.0800 0 0 4.0319 2.6879 1.6800 0.9600 0.4800 0.1920 0.0480 0 0 0 4.8383 2.6880 1.3440 0.5760 0.1920 0.0384 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cond(Ahat)

ans = 7.2583e+04

[Qiang Li]

Hi John,

The full matrices are

K>> num2str(A_2stepblock,100)

ans =

  10×1003 char array

    '                                       0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          1'
    '                                       0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          1                                                                                                          0'
    '                                       0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          2                                                                                                          0                                                                                                          0'
    '                                       0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          6                                                                                                          0                                                                                                          0                                                                                                          0'
    '                                       0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                         24                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'
    '1953050.78182335919700562953948974609375                                                                   390611.805629631155170500278472900390625                                                                   78122.6909803108428604900836944580078125                                                              15624.604167217634312692098319530487060546875                                                              3124.9340277303299444611184298992156982421875                                                             624.989444414568652064190246164798736572265625                                                         124.9984166588426575117409811355173587799072265625                                                           24.999788887400001158312079496681690216064453125                                                       4.99997888869543150036633960553444921970367431640625                                                                                                          1'
    ' 3515506.2506666802801191806793212890625                                                                     624981.5278424867428839206695556640625                                                                 109372.22917052343836985528469085693359375                                                                 18749.604166381977847777307033538818359375                                                              3124.9472220728430329472757875919342041015625                                                           499.99366663537063004696392454206943511962890625                                                           74.999366662200003474936238490045070648193359375                                                        9.9999577773908630007326792110688984394073486328125                                                                                                          1                                                                                                          0'
    '   5624833.75058238022029399871826171875                                                                    874977.83336418750695884227752685546875                                                                 131247.22916467385948635637760162353515625                                                                18749.6833324370600166730582714080810546875                                                             2499.96833317685332076507620513439178466796875                                                            299.9974666488000138997449539601802825927734375                                                         29.99987333217259077855487703345715999603271484375                                                                                                          2                                                                                                          0                                                                                                          0'
    '  7874800.500277687795460224151611328125                                                                   1049977.83331739087589085102081298828125                                                                  131247.7833270594128407537937164306640625                                                               14999.80999906111901509575545787811279296875                                                                1499.98733324400018318556249141693115234375                                                          119.999493328690363114219508133828639984130859375                                                                                                          6                                                                                                          0                                                                                                          0                                                                                                          0'
    '   9449800.49985651858150959014892578125                                                                    1049982.2666164753027260303497314453125                                                                104998.669993427829467691481113433837890625                                                                 8999.9239994640010991133749485015869140625                                                            599.9974666434517303059692494571208953857421875                                                                                                         24                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'

K>> num2str(dA_2stepblock_dt,100)

ans =

  10×1002 char array

    '                                      0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'
    '                                      0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'
    '                                      0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'
    '                                      0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'
    '                                      0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'
    '3515506.2506666802801191806793212890625                                                                     624981.5278424867428839206695556640625                                                                 109372.22917052343836985528469085693359375                                                                 18749.604166381977847777307033538818359375                                                              3124.9472220728430329472757875919342041015625                                                           499.99366663537063004696392454206943511962890625                                                           74.999366662200003474936238490045070648193359375                                                        9.9999577773908630007326792110688984394073486328125                                                                                                          1                                                                                                          0'
    '  5624833.75058238022029399871826171875                                                                    874977.83336418750695884227752685546875                                                                 131247.22916467385948635637760162353515625                                                                18749.6833324370600166730582714080810546875                                                             2499.96833317685332076507620513439178466796875                                                            299.9974666488000138997449539601802825927734375                                                         29.99987333217259077855487703345715999603271484375                                                                                                          2                                                                                                          0                                                                                                          0'
    ' 7874800.500277687795460224151611328125                                                                   1049977.83331739087589085102081298828125                                                                  131247.7833270594128407537937164306640625                                                               14999.80999906111901509575545787811279296875                                                                1499.98733324400018318556249141693115234375                                                          119.999493328690363114219508133828639984130859375                                                                                                          6                                                                                                          0                                                                                                          0                                                                                                          0'
    '  9449800.49985651858150959014892578125                                                                    1049982.2666164753027260303497314453125                                                                104998.669993427829467691481113433837890625                                                                 8999.9239994640010991133749485015869140625                                                            599.9974666434517303059692494571208953857421875                                                                                                         24                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'
    '     9449840.39954827725887298583984375                                                                   839989.359947422635741531848907470703125                                                                 62999.467996248000417836010456085205078125                                                             3599.98479986071060920949093997478485107421875                                                                                                        120                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0                                                                                                          0'

I honestly don't know how to copy this out properly, but hope this will give you a better idea.

[Qiang Li]

A = A_2stepblock

B = dA_2stepblock_dt

[Torsten]

I honestly don't know how to copy this out properly, but hope this will give you a better idea.

Save them in a .mat file and include this file here.

You should stop working from the command line. Use scripts and functions.

[Qiang Li]

This is how they are constructed

time_points =[ 0
   4.999978888695432
   9.999771309046501
  14.999222028535593
  19.999134717433272
  24.999047406330948
  27.468441858944786];

A_2stepblock = [A(time_points(i_seg));A(time_points(i_seg+1))];

function output = A(time)
%build the A matrix at a time point
output = [     time^9,      time^8,     time^7,     time^6,    time^5,    time^4,   time^3, time^2, time, 1
   9*time^8,    8*time^7,   7*time^6,   6*time^5,  5*time^4,  4*time^3, 3*time^2, 2*time,       1, 0
  72*time^7,   56*time^6,  42*time^5,  30*time^4, 20*time^3, 12*time^2,   6*time,         2,       0, 0
 504*time^6,  336*time^5, 210*time^4, 120*time^3, 60*time^2,   24*time,           6,         0,       0, 0
3024*time^5, 1680*time^4, 840*time^3, 360*time^2,  120*time,           24,           0,         0,       0, 0];
end

and

dA_2stepblock_dt = [zeros(10,5);dAdt(time_points(i_seg+1))];

function output = dAdt(time)
%build the dAdt matrix at a time point
output = [   9*time^8,    8*time^7,   7*time^6,   6*time^5,  5*time^4,  4*time^3, 3*time^2, 2*time,       1, 0
  72*time^7,   56*time^6,  42*time^5,  30*time^4, 20*time^3, 12*time^2,   6*time,         2,       0, 0
 504*time^6,  336*time^5, 210*time^4, 120*time^3, 60*time^2,   24*time,           6,         0,       0, 0
3024*time^5, 1680*time^4, 840*time^3, 360*time^2,  120*time,           24,           0,         0,       0, 0
15120*time^4,   6720*time^3,    2520*time^2,    720*time,   120,    0,  0,  0,  0,  0];
end

For matrices in the question, I used i_seg = 1.

[Qiang Li]

I'm not familiar with decomposition and re-scaling as you mentioned. Can you provide more on how they could help here? Thanks.

[John D'Errico]

Decomposition effectively computes the LU factorization of your matrix A, storing it into the object I called Afac. But then instead of performing solves, it uses back and forward substitutions, since L and U lower and upper are triangular matrices respectively. I have a funny feeling that the direct slash calls will be pretty much identical in their result as LU, since those tools will also generate an internal decomposition of the same form, and using the explicit decomposition tool will probably have more overhead.

As far as scaling goes, does it matter if you measure time in seconds, or in milliseconds? NO!!!!!! Of course not. However, it DOES matter when you perform linear algebra. For example, a simple Vandermonde matrix, used to perform things like polynomial estimation, interpolation, etc., is a great example.

x = sort(rand(10,1)*27)

x = 10x1

7.4456 8.6356 13.7367 17.9532 20.4170 20.4825 20.6787 23.3689 24.8598 26.5335

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Now, build a vandermonde matrix.

V = x.^(9:-1:0);
cond(V)

ans = 4.6779e+20

Now, just scale the vector x. This is no different than a different choice of units. Does it matter if I work in terms of meters, versus centimeters?

xhat = x/25

xhat = 10x1

0.2978 0.3454 0.5495 0.7181 0.8167 0.8193 0.8271 0.9348 0.9944 1.0613

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

Vhat = xhat.^(9:-1:0);
cond(Vhat)

ans = 6.5979e+11

All it took to reduce the condition number of that matrix by a factor of 1e8, is to rescale the data, effectively, changing the units on my vector x.

Now, look at your problem.

time_points =[ 0
    4.999978888695432
    9.999771309046501
    14.999222028535593
    19.999134717433272
    24.999047406330948
    27.468441858944786];

time_scaled = time_points/5;

i_seg = 1;

A_2stepblock = [A(time_scaled(i_seg));A(time_scaled(i_seg+1))]

A_2stepblock = 10x10

1.0e+03 * 0 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0.0010 0 0 0 0 0 0 0 0 0.0020 0 0 0 0 0 0 0 0 0.0060 0 0 0 0 0 0 0 0 0.0240 0 0 0 0 0.0010 0.0010 0.0010 0.0010 0.0010 0.0010 0.0010 0.0010 0.0010 0.0010 0.0090 0.0080 0.0070 0.0060 0.0050 0.0040 0.0030 0.0020 0.0010 0 0.0720 0.0560 0.0420 0.0300 0.0200 0.0120 0.0060 0.0020 0 0 0.5040 0.3360 0.2100 0.1200 0.0600 0.0240 0.0060 0 0 0 3.0239 1.6800 0.8400 0.3600 0.1200 0.0240 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

cond(A_2stepblock)

ans = 4.4345e+06

Disappointingly, it does not massively decrease the condition number, but even that extent of decrease may be significant, and sufficient to help.

function output = A(time)
%build the A matrix at a time point
output = [     time^9,      time^8,     time^7,     time^6,    time^5,    time^4,   time^3, time^2, time, 1
    9*time^8,    8*time^7,   7*time^6,   6*time^5,  5*time^4,  4*time^3, 3*time^2, 2*time,       1, 0
    72*time^7,   56*time^6,  42*time^5,  30*time^4, 20*time^3, 12*time^2,   6*time,         2,       0, 0
    504*time^6,  336*time^5, 210*time^4, 120*time^3, 60*time^2,   24*time,           6,         0,       0, 0
    3024*time^5, 1680*time^4, 840*time^3, 360*time^2,  120*time,           24,           0,         0,       0, 0];
end

% dA_2stepblock_dt = [zeros(10,5);dAdt(time_points(i_seg+1))];

function output = dAdt(time)
%build the dAdt matrix at a time point
output = [   9*time^8,    8*time^7,   7*time^6,   6*time^5,  5*time^4,  4*time^3, 3*time^2, 2*time,       1, 0
    72*time^7,   56*time^6,  42*time^5,  30*time^4, 20*time^3, 12*time^2,   6*time,         2,       0, 0
    504*time^6,  336*time^5, 210*time^4, 120*time^3, 60*time^2,   24*time,           6,         0,       0, 0
    3024*time^5, 1680*time^4, 840*time^3, 360*time^2,  120*time,           24,           0,         0,       0, 0
    15120*time^4,   6720*time^3,    2520*time^2,    720*time,   120,    0,  0,  0,  0,  0];
end

[Qiang Li]

John, in your original answer, the diag(5.^(-9:1:0)) rescaling worked really well.

Ahat = [A(1:5,:);A(6:10,:)*diag(5.^(-9:1:0))]

Can you tell me how you choose this rescaling matrix? How am I supposed to transform the scale back after inversing Ahat in this case? Thanks.

[Qiang Li]

I think for small matrices as in my case, decomposing the matrix is not worth it, as decomposition costs significantly more than the savings from the inverse. It may work for large matrices though.

[Qiang Li]

The rescaling suggestion is making a lot of sense now. Definitely need to rescale A to improve the accuracy.

Answer 2

[Matt J]

I tried both inv(A)*B*inv(A) or A\B/A, resulting in very different results for a gradient descent algorithm.

From this, it would appear that what you are really trying to do is solve A*inv(B)*A*x=y. If so, you should consider lscov,

x=lscov(A,y,B);

[Qiang Li]

Hi Matt,

I'm not solving for x at the moment. But I may actually need to do it later in my project. So thanks in advance!

Answer 3

[Qiang Li]

The short answer for my application: A\B/A.

inv(A)*B*inv(A) introduces inconsistency after iterations.

[Qiang Li]

Neither decomposition or exploiting the structure of matrix A improves the inconsistency.