Thread Subject: inverse matrix inv pinv linfactor ginv

  
Hi there,

I want to solve or b the : (X'X)*b=X'*y. I have some issues with the inverse of X'*X.
(vector b should estimate the vector [5 0.1048]. (regression coefficients))
MATLAB(R2007b).

In the beginning I used "inv" and got a warning:

Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 1.000436e-148.

b =

    4.9346
    0.0000 (Inspitethe warning the solution seems to approach the real values)

I also used pinv but I am not sure about the result
b =

  1.0e-073 *

    0.0000
    0.8577 (Here the result looks weird, don't get a warning though)

I also used linfactor (Timothy A. Davis)

bnaive=linfactor (linfactor(Ximp'*Ximp),Ximp'*y')
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 3.483766e-074.
> In linfactor at 175
Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 3.483766e-074.
> In linfactor at 175

b =

    4.9346
    0.0000 (Again a warning but a good looking solution)

And finally ginv (Luis Frank, jan 2009.)

and got either [5.4417 0] or [1.0e-072 * 0 0.1171] without warnings.

So I am a little confused..Which is correct?...Below I give the data. Note that X=[ones(n,1) x]. where n=300.


x=
1.0e+073 *

    6.0217
    6.9355
    6.9355
    0.0000
    6.9355
    6.9355
    0.0000
    6.6810
    4.5431
    6.1384
    0.0000
    5.8016
    6.2875
    6.9355
    4.5401
    0.0000
    6.9355
    6.6662
    6.9355
    5.9197
    6.9355
    0.0000
    6.9355
    6.1304
    5.6790
    6.9355
    6.9355
    6.9355
    6.9355
    4.7936
    0.0000
    0.0000
    6.6226
    0.0000
    0.0000
    0.0000
    6.9355
    6.9355
    4.5560
    4.9384
    6.9355
    5.7786
    0.0000
    6.9355
    0.0000
    5.4005
    6.2339
    6.9355
    0.0000
    0.0000
    5.1291
    6.9355
    6.9355
    6.9355
    6.9355
    0.0000
    0.0000
    6.9355
    6.9355
    0.0000
    5.1252
    6.1170
    0.0000
    6.9355
    6.9355
    6.9355
    6.9355
    0.0000
    5.0401
    6.8300
    4.8740
    0.0000
    6.6515
    0.0000
    5.0013
    0.0000
    6.9355
    0.0000
    4.5467
    5.7678
    6.9355
    0.0000
    6.5277
    4.9925
    6.9355
    6.1809
    5.0394
    5.2346
    0.0000
    0.0000
    6.8909
    6.9355
    5.6816
    4.9700
    0.0000
    0.0000
    5.0839
    4.9574
    0.0000
    5.2523
    5.1478
    6.9355
    6.9355
    0.0000
    0.0000
    4.5564
    5.8503
    0.0000
    6.9355
    4.5315
    4.6217
    5.7917
    6.9355
    0.0000
    6.9355
    5.0533
    6.9355
    6.9355
    0.0000
    6.9355
    0.0000
    6.9355
    6.3646
    6.9355
    0.0000
    5.4656
    0.0000
    6.9355
    6.9355
    6.9355
    6.9355
    0.0000
    0.0000
    4.8783
    6.9355
    6.9355
    0.0000
    4.7015
    4.8418
    4.8581
    6.9355
    6.9355
    6.9355
    5.4954
    0.0000
    6.3505
    6.9355
    6.9355
    6.9355
    0.0000
    6.9355
    0.0000
    6.9355
    0.0000
    5.9320
    6.9355
    4.8874
    6.9355
    6.4735
    6.9355
    6.9355
    0.0000
    5.6360
    5.2314
    6.9355
    6.9355
    5.5291
    6.9355
    6.8946
    6.9355
    6.9355
    0.0000
    6.7673
    5.3281
    0.0000
    6.9355
    6.9355
    0.0000
    6.9355
    5.9098
    6.9355
    6.9355
    0.0000
    0.0000
    6.9355
    6.9355
    0.0000
    6.9355
    0.0000
    5.3169
    5.4155
    0.0000
    6.9355
    6.0202
    0.0000
    0.0000
    6.9355
    6.9355
    6.9355
    5.7784
    5.9336
    6.9355
    5.6023
    0.0000
    6.9355
    0.0000
    6.9355
    4.5746
    6.9355
    0.0000
    6.9355
    4.5628
    0.0000
    0.0000
    4.6689
    5.9038
    0.0000
    4.9554
    4.9408
    6.9355
    4.5740
    4.7817
    6.9355
    5.6149
    6.9355
    4.7791
    4.5862
    5.8192
    0.0000
    0.0000
    4.5331
    5.2861
    6.9355
    6.9355
    5.4413
    6.9355
    5.3689
    0.0000
    4.6382
    5.2439
    0.0000
    6.9355
    0.0000
    6.9355
    6.9355
    5.2391
    6.9355
    6.9355
    6.9355
    0.0000
    6.9355
    6.9355
    5.6362
    6.9355
    0.0000
    0.0000
    0.0000
    6.9355
    6.9355
    6.8841
    0.0000
    5.0142
    6.9355
    5.4128
    0.0000
    0.0000
    0.0000
    6.0928
    6.5998
    6.9355
    0.0000
    6.9355
    4.7202
    6.9355
    4.7974
    5.1762
    0.0000
    6.9355
    6.5298
    0.0000
    0.0000
    6.9355
    6.1336
    5.4558
    0.0000
    0.0000
    5.6442
    0.0000
    5.2961
    6.9355
    5.2238
    6.6430
    6.9355
    6.9355
    6.9355
    0.0000
    6.9355
    5.2214
    6.9355
    5.3165





y=

    4.8111
    5.2604
    6.3000
    4.1397
    5.7850
    8.0680
    2.6153
    7.1341
    6.3813
    5.1344
    4.8746
    6.0043
    5.6167
    5.6686
    5.7879
    4.9851
    3.6736
    5.5030
    5.2088
    4.8137
    4.5450
    5.3842
    6.1297
    5.1319
    4.4700
    4.7678
    6.2697
    5.3731
    5.0480
    6.1615
    5.4038
    4.4000
    4.7029
    6.3949
    5.6621
    4.0477
    3.5602
    4.6434
    3.3360
    6.5198
    4.1286
    6.4947
    4.7885
    7.4672
    5.7104
    6.4073
    5.3284
    5.6730
    4.6356
    4.8480
    4.2965
    5.5603
    4.5348
    3.5832
    5.4171
    3.9930
    4.4879
    4.5334
    6.2194
    5.2592
    3.0510
    5.2886
    5.4761
    5.2904
    7.3873
    5.0014
    4.8700
    6.0900
    6.1828
    4.2469
    6.5676
    6.3699
    6.2629
    4.0863
    6.6842
    5.3526
    5.6865
    4.5678
    5.3424
    5.1886
    4.5782
    4.1196
    4.4998
    3.3030
    4.5550
    6.1645
    4.8040
    6.3907
    4.4269
    4.4374
    5.4680
    5.7414
    7.3257
    6.1399
    2.9633
    4.5267
    5.0329
    5.7673
    4.0117
    5.6086
    5.6665
    4.2401
    4.2609
    5.3666
    4.9009
    5.3461
    6.6157
    5.3064
    3.0179
    5.4734
    5.8616
    2.4153
    7.7829
    3.8170
    5.0442
    5.2099
    5.3179
    4.8733
    4.2587
    6.2033
    5.1428
    5.3150
    7.5840
    5.6033
    5.5261
    6.4733
    6.0306
    4.9166
    4.4991
    4.8719
    5.5547
    6.5007
    5.8086
    4.2790
    5.7095
    5.4276
    6.6512
    5.2374
    7.2476
    4.2253
    4.7902
    7.4058
    6.3936
    5.5561
    5.6328
    7.4945
    7.6635
    5.7805
    5.1112
    2.7504
    4.1025
    6.9345
    5.1199
    4.9941
    4.5647
    5.4471
    4.4141
    6.7696
    4.7769
    5.2200
    6.0108
    3.3302
    7.3792
    5.5906
    4.9745
    5.3593
    5.9843
    5.7817
    5.0402
    6.0475
    5.5631
    5.7827
    4.4960
    5.1135
    3.3456
    4.6051
    6.6434
    5.9896
    5.9918
    5.3772
    5.5311
    4.5358
    5.0157
    4.4853
    6.6768
    7.7551
    5.4937
    6.8068
    4.5231
    4.2895
    5.0643
    4.9405
    4.0419
    5.7608
    5.1294
    5.0549
    6.1442
    3.8324
    5.8736
    6.0473
    6.0278
    5.8280
    4.7676
    5.5682
    3.9522
    4.2703
    5.3895
    8.2634
    4.2123
    4.6494
    5.7650
    5.0053
    5.8436
    4.5409
    6.0050
    7.0507
    5.7757
    5.6222
    5.6440
    5.8824
    2.8991
    4.4518
    3.7672
    5.0790
    4.9698
    4.8897
    5.7037
    6.3800
    5.2329
    5.6111
    4.8948
    6.2990
    3.1272
    4.5467
    4.2359
    4.1112
    5.6829
    4.6515
    5.6744
    4.5980
    3.9540
    5.2330
    4.4072
    6.5943
    4.3242
    4.0208
    6.2833
    6.7559
    5.9868
    4.6251
    7.6454
    6.0673
    5.8195
    6.1676
    4.6909
    2.9627
    5.2211
    5.4769
    5.8216
    4.6540
    6.1510
    5.9472
    6.8202
    5.9406
    7.1331
    4.8461
    4.9920
    7.4089
    5.8690
    6.6683
    3.4887
    3.5328
    4.7221
    5.7539
    4.4846
    5.2311
    6.5396
    6.0771
    5.4340
    4.1203
    5.8677
    6.3185
    3.8214
    6.7323
    4.0473
    2.8361
    5.0833
    4.1271
    5.7618
    3.8940
    4.7097
    4.9934
    3.7886
    4.9597
    6.5558
    5.7036
    4.1438
    5.5269
    6.2652
    7.1042

"leo nidas" <bleonidas25@yahoo.gr> wrote in message <ha8r6f$pf7$1@fred.mathworks.com>...
>
> Hi there,
>
> I want to solve or b the : (X'X)*b=X'*y. I have some issues with the inverse of X'*X.
> (vector b should estimate the vector [5 0.1048]. (regression coefficients))
> MATLAB(R2007b).
>
> In the beginning I used "inv" and got a warning:
>
> Warning: Matrix is close to singular or badly scaled.
> Results may be inaccurate. RCOND = 1.000436e-148.
>
> b =
>
> 4.9346
> 0.0000 (Inspitethe warning the solution seems to approach the real values)
>
> I also used pinv but I am not sure about the result
> b =
>
> 1.0e-073 *
>
> 0.0000
> 0.8577 (Here the result looks weird, don't get a warning though)
>
> I also used linfactor (Timothy A. Davis)
>
> bnaive=linfactor (linfactor(Ximp'*Ximp),Ximp'*y')
> Warning: Matrix is close to singular or badly scaled.
> Results may be inaccurate. RCOND = 3.483766e-074.
> > In linfactor at 175
> Warning: Matrix is close to singular or badly scaled.
> Results may be inaccurate. RCOND = 3.483766e-074.
> > In linfactor at 175
>
> b =
>
> 4.9346
> 0.0000 (Again a warning but a good looking solution)
>
> And finally ginv (Luis Frank, jan 2009.)
>
> and got either [5.4417 0] or [1.0e-072 * 0 0.1171] without warnings.
>
> So I am a little confused..Which is correct?...Below I give the data. Note that X=[ones(n,1) x]. where n=300.
>
>
> x=
> 1.0e+073 *
>
> 6.0217
> 6.9355
> 6.9355
> 0.0000
> 6.9355
> 6.9355
> 0.0000
> 6.6810
> 4.5431
> 6.1384
> 0.0000
> 5.8016
> 6.2875
> 6.9355
> 4.5401
> 0.0000
> 6.9355
> 6.6662
> 6.9355
> 5.9197
> 6.9355
> 0.0000
> 6.9355
> 6.1304
> 5.6790
> 6.9355
> 6.9355
> 6.9355
> 6.9355
> 4.7936
> 0.0000
> 0.0000
> 6.6226
> 0.0000
> 0.0000
> 0.0000
> 6.9355
> 6.9355
> 4.5560
> 4.9384
> 6.9355
> 5.7786
> 0.0000
> 6.9355
> 0.0000
> 5.4005
> 6.2339
> 6.9355
> 0.0000
> 0.0000
> 5.1291
> 6.9355
> 6.9355
> 6.9355
> 6.9355
> 0.0000
> 0.0000
> 6.9355
> 6.9355
> 0.0000
> 5.1252
> 6.1170
> 0.0000
> 6.9355
> 6.9355
> 6.9355
> 6.9355
> 0.0000
> 5.0401
> 6.8300
> 4.8740
> 0.0000
> 6.6515
> 0.0000
> 5.0013
> 0.0000
> 6.9355
> 0.0000
> 4.5467
> 5.7678
> 6.9355
> 0.0000
> 6.5277
> 4.9925
> 6.9355
> 6.1809
> 5.0394
> 5.2346
> 0.0000
> 0.0000
> 6.8909
> 6.9355
> 5.6816
> 4.9700
> 0.0000
> 0.0000
> 5.0839
> 4.9574
> 0.0000
> 5.2523
> 5.1478
> 6.9355
> 6.9355
> 0.0000
> 0.0000
> 4.5564
> 5.8503
> 0.0000
> 6.9355
> 4.5315
> 4.6217
> 5.7917
> 6.9355
> 0.0000
> 6.9355
> 5.0533
> 6.9355
> 6.9355
> 0.0000
> 6.9355
> 0.0000
> 6.9355
> 6.3646
> 6.9355
> 0.0000
> 5.4656
> 0.0000
> 6.9355
> 6.9355
> 6.9355
> 6.9355
> 0.0000
> 0.0000
> 4.8783
> 6.9355
> 6.9355
> 0.0000
> 4.7015
> 4.8418
> 4.8581
> 6.9355
> 6.9355
> 6.9355
> 5.4954
> 0.0000
> 6.3505
> 6.9355
> 6.9355
> 6.9355
> 0.0000
> 6.9355
> 0.0000
> 6.9355
> 0.0000
> 5.9320
> 6.9355
> 4.8874
> 6.9355
> 6.4735
> 6.9355
> 6.9355
> 0.0000
> 5.6360
> 5.2314
> 6.9355
> 6.9355
> 5.5291
> 6.9355
> 6.8946
> 6.9355
> 6.9355
> 0.0000
> 6.7673
> 5.3281
> 0.0000
> 6.9355
> 6.9355
> 0.0000
> 6.9355
> 5.9098
> 6.9355
> 6.9355
> 0.0000
> 0.0000
> 6.9355
> 6.9355
> 0.0000
> 6.9355
> 0.0000
> 5.3169
> 5.4155
> 0.0000
> 6.9355
> 6.0202
> 0.0000
> 0.0000
> 6.9355
> 6.9355
> 6.9355
> 5.7784
> 5.9336
> 6.9355
> 5.6023
> 0.0000
> 6.9355
> 0.0000
> 6.9355
> 4.5746
> 6.9355
> 0.0000
> 6.9355
> 4.5628
> 0.0000
> 0.0000
> 4.6689
> 5.9038
> 0.0000
> 4.9554
> 4.9408
> 6.9355
> 4.5740
> 4.7817
> 6.9355
> 5.6149
> 6.9355
> 4.7791
> 4.5862
> 5.8192
> 0.0000
> 0.0000
> 4.5331
> 5.2861
> 6.9355
> 6.9355
> 5.4413
> 6.9355
> 5.3689
> 0.0000
> 4.6382
> 5.2439
> 0.0000
> 6.9355
> 0.0000
> 6.9355
> 6.9355
> 5.2391
> 6.9355
> 6.9355
> 6.9355
> 0.0000
> 6.9355
> 6.9355
> 5.6362
> 6.9355
> 0.0000
> 0.0000
> 0.0000
> 6.9355
> 6.9355
> 6.8841
> 0.0000
> 5.0142
> 6.9355
> 5.4128
> 0.0000
> 0.0000
> 0.0000
> 6.0928
> 6.5998
> 6.9355
> 0.0000
> 6.9355
> 4.7202
> 6.9355
> 4.7974
> 5.1762
> 0.0000
> 6.9355
> 6.5298
> 0.0000
> 0.0000
> 6.9355
> 6.1336
> 5.4558
> 0.0000
> 0.0000
> 5.6442
> 0.0000
> 5.2961
> 6.9355
> 5.2238
> 6.6430
> 6.9355
> 6.9355
> 6.9355
> 0.0000
> 6.9355
> 5.2214
> 6.9355
> 5.3165
>
>
>
>
>
> y=
>
> 4.8111
> 5.2604
> 6.3000
> 4.1397
> 5.7850
> 8.0680
> 2.6153
> 7.1341
> 6.3813
> 5.1344
> 4.8746
> 6.0043
> 5.6167
> 5.6686
> 5.7879
> 4.9851
> 3.6736
> 5.5030
> 5.2088
> 4.8137
> 4.5450
> 5.3842
> 6.1297
> 5.1319
> 4.4700
> 4.7678
> 6.2697
> 5.3731
> 5.0480
> 6.1615
> 5.4038
> 4.4000
> 4.7029
> 6.3949
> 5.6621
> 4.0477
> 3.5602
> 4.6434
> 3.3360
> 6.5198
> 4.1286
> 6.4947
> 4.7885
> 7.4672
> 5.7104
> 6.4073
> 5.3284
> 5.6730
> 4.6356
> 4.8480
> 4.2965
> 5.5603
> 4.5348
> 3.5832
> 5.4171
> 3.9930
> 4.4879
> 4.5334
> 6.2194
> 5.2592
> 3.0510
> 5.2886
> 5.4761
> 5.2904
> 7.3873
> 5.0014
> 4.8700
> 6.0900
> 6.1828
> 4.2469
> 6.5676
> 6.3699
> 6.2629
> 4.0863
> 6.6842
> 5.3526
> 5.6865
> 4.5678
> 5.3424
> 5.1886
> 4.5782
> 4.1196
> 4.4998
> 3.3030
> 4.5550
> 6.1645
> 4.8040
> 6.3907
> 4.4269
> 4.4374
> 5.4680
> 5.7414
> 7.3257
> 6.1399
> 2.9633
> 4.5267
> 5.0329
> 5.7673
> 4.0117
> 5.6086
> 5.6665
> 4.2401
> 4.2609
> 5.3666
> 4.9009
> 5.3461
> 6.6157
> 5.3064
> 3.0179
> 5.4734
> 5.8616
> 2.4153
> 7.7829
> 3.8170
> 5.0442
> 5.2099
> 5.3179
> 4.8733
> 4.2587
> 6.2033
> 5.1428
> 5.3150
> 7.5840
> 5.6033
> 5.5261
> 6.4733
> 6.0306
> 4.9166
> 4.4991
> 4.8719
> 5.5547
> 6.5007
> 5.8086
> 4.2790
> 5.7095
> 5.4276
> 6.6512
> 5.2374
> 7.2476
> 4.2253
> 4.7902
> 7.4058
> 6.3936
> 5.5561
> 5.6328
> 7.4945
> 7.6635
> 5.7805
> 5.1112
> 2.7504
> 4.1025
> 6.9345
> 5.1199
> 4.9941
> 4.5647
> 5.4471
> 4.4141
> 6.7696
> 4.7769
> 5.2200
> 6.0108
> 3.3302
> 7.3792
> 5.5906
> 4.9745
> 5.3593
> 5.9843
> 5.7817
> 5.0402
> 6.0475
> 5.5631
> 5.7827
> 4.4960
> 5.1135
> 3.3456
> 4.6051
> 6.6434
> 5.9896
> 5.9918
> 5.3772
> 5.5311
> 4.5358
> 5.0157
> 4.4853
> 6.6768
> 7.7551
> 5.4937
> 6.8068
> 4.5231
> 4.2895
> 5.0643
> 4.9405
> 4.0419
> 5.7608
> 5.1294
> 5.0549
> 6.1442
> 3.8324
> 5.8736
> 6.0473
> 6.0278
> 5.8280
> 4.7676
> 5.5682
> 3.9522
> 4.2703
> 5.3895
> 8.2634
> 4.2123
> 4.6494
> 5.7650
> 5.0053
> 5.8436
> 4.5409
> 6.0050
> 7.0507
> 5.7757
> 5.6222
> 5.6440
> 5.8824
> 2.8991
> 4.4518
> 3.7672
> 5.0790
> 4.9698
> 4.8897
> 5.7037
> 6.3800
> 5.2329
> 5.6111
> 4.8948
> 6.2990
> 3.1272
> 4.5467
> 4.2359
> 4.1112
> 5.6829
> 4.6515
> 5.6744
> 4.5980
> 3.9540
> 5.2330
> 4.4072
> 6.5943
> 4.3242
> 4.0208
> 6.2833
> 6.7559
> 5.9868
> 4.6251
> 7.6454
> 6.0673
> 5.8195
> 6.1676
> 4.6909
> 2.9627
> 5.2211
> 5.4769
> 5.8216
> 4.6540
> 6.1510
> 5.9472
> 6.8202
> 5.9406
> 7.1331
> 4.8461
> 4.9920
> 7.4089
> 5.8690
> 6.6683
> 3.4887
> 3.5328
> 4.7221
> 5.7539
> 4.4846
> 5.2311
> 6.5396
> 6.0771
> 5.4340
> 4.1203
> 5.8677
> 6.3185
> 3.8214
> 6.7323
> 4.0473
> 2.8361
> 5.0833
> 4.1271
> 5.7618
> 3.8940
> 4.7097
> 4.9934
> 3.7886
> 4.9597
> 6.5558
> 5.7036
> 4.1438
> 5.5269
> 6.2652
> 7.1042


use the backslash operator:
b=(X'X)\X'y;

see ex. here:
http://blogs.mathworks.com/pick/2009/06/26/dont-let-that-inv-go-past-your-eyes-to-solve-that-system-factorize/

and here: http://blogs.mathworks.com/loren/2007/05/16/purpose-of-inv/

I forgot to mention that I already have tried the backslash operator and I get a warning :

Warning: Matrix is close to singular or badly scaled.
         Results may be inaccurate. RCOND = 1.000436e-148.

ans =

    4.9346
    0.0000

Thanx for the answer! And for any other answers in advance.

"leo nidas" <bleonidas25@yahoo.gr> wrote in message <ha8r6f$pf7$1@fred.mathworks.com>...


 Note that X=[ones(n,1) x]. where n=300.
>

When you do linear algebra, the first thing is to SCALE your variable correctly. Combine a vector of a magnitude 1 with x in the order of 1e73 is very bad.

I assume you want to find a linear polynomial coefficients that fit your data.

In this case work with xi := x/1e73 (xi will be between 0 and 1)

The fitting you want
y ~= b0 + b1*x becomes
y ~= [ones(n,1) x/1e73] * [b0 b1*1e73]'

solve for b2 = [b0 b1*1e73] by least square

xi = x/1e73
A2 = [ones(n,1) xi]
b2 = A2 \ y
b = [b2(1) b2(2)/1e73]

% Bruno

PS: A quick look of GINV indicates it's a not good quality code, I would say even amateurish one.

"Petros Mpogiatzis" <painter@geo.auth.gr> wrote in message <ha8s85$4qh$1@fred.mathworks.com>...
> "leo nidas" <bleonidas25@yahoo.gr> wrote in message <ha8r6f$pf7$1@fred.mathworks.com>...
> >
> > Hi there,
> >
> > I want to solve or b the : (X'X)*b=X'*y. I have some issues with the inverse of X'*X.
> > (vector b should estimate the vector [5 0.1048]. (regression coefficients))
> > MATLAB(R2007b).
> >
> > In the beginning I used "inv" and got a warning:
> >
> > Warning: Matrix is close to singular or badly scaled.
> > Results may be inaccurate. RCOND = 1.000436e-148.
> >
> > b =
> >
> > 4.9346
> > 0.0000 (Inspitethe warning the solution seems to approach the real values)
> >
> > I also used pinv but I am not sure about the result
> > b =
> >
> > 1.0e-073 *
> >
> > 0.0000
> > 0.8577 (Here the result looks weird, don't get a warning though)
> >
> > I also used linfactor (Timothy A. Davis)
> >
> > bnaive=linfactor (linfactor(Ximp'*Ximp),Ximp'*y')
> > Warning: Matrix is close to singular or badly scaled.
> > Results may be inaccurate. RCOND = 3.483766e-074.
> > > In linfactor at 175
> > Warning: Matrix is close to singular or badly scaled.
> > Results may be inaccurate. RCOND = 3.483766e-074.
> > > In linfactor at 175
> >
> > b =
> >
> > 4.9346
> > 0.0000 (Again a warning but a good looking solution)
> >
> > And finally ginv (Luis Frank, jan 2009.)
> >
> > and got either [5.4417 0] or [1.0e-072 * 0 0.1171] without warnings.
> >
> > So I am a little confused..Which is correct?...Below I give the data. Note that X=[ones(n,1) x]. where n=300.
> >
> >
> > x=
> > 1.0e+073 *
> >
> > 6.0217

(snip)

> > 6.2652
> > 7.1042
>
>
> use the backslash operator:
> b=(X'X)\X'y;

NO!!!!!

Don't just blindly use backslash on the normal equations!!!!!
Once you form the normal equations, you have already
screwed yourself.

Sigh.

Don't use the normal equations to solve a linear regression
problem in the first place. Solve it as

  b = X\y;

Or, solve it as

  b = pinv(X)*y;

This is the correct way to solve the linear regression problem.
When you use the normal equations, as people seem to teach
each other, INCORRECTLY, to solve a linear regression problem,
you worsen the conditioning of any regression that you will
solve.

But BEFORE you do that, scale this data properly, as I show
how below.

So next, we must deal with that little, inconsequential factor
of 1e73 on your data! Look at the variable x! See the 1e73
that comes out front. It says

> > x=
> > 1.0e+073 *
> >
> > 6.0217

This means that this regression will generally return garbage,
unless you are careful.

Don't just throw such unscaled data into a linear algebra
problem, and worse, then throw the normal equations at it!
It is this that is causing the other part of your problem.

b = [ones(300,1),x/1e73]\y;
b =
       4.9346
     0.079934

Only now should you rescale b to compensate for the
prior scaling of x.

format long g
b(2) = b(2) * 1e73
b =
          4.93457076016252
      7.99340988034286e+71

See that no warning message was generated now.

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <ha9rl7$gej$1@fred.mathworks.com>...

>
> Only now should you rescale b to compensate for the
> prior scaling of x.
>
> format long g
> b(2) = b(2) * 1e73
> b =
> 4.93457076016252


Of course, John meant divided b(2) is a correct rescaling.

Bruno

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <ha9sh9$l13$1@fred.mathworks.com>...
> "John D'Errico" <woodchips@rochester.rr.com> wrote in message <ha9rl7$gej$1@fred.mathworks.com>...
>
> >
> > Only now should you rescale b to compensate for the
> > prior scaling of x.
> >
> > format long g
> > b(2) = b(2) * 1e73
> > b =
> > 4.93457076016252
>
>
> Of course, John meant divided b(2) is a correct rescaling.
>
> Bruno

What is a little factor of 1e73 one way or the other among
friends anyway?

John

See my FACTORIZE code; it's very easy to use. Also, posting a matrix with 5 digits of accuracy isn't a good idea. It's too long for anyone to read, and even if someone wants to reconstruct your problem, they can't because you've tossed out about 11 digits in each number (which is where the problem may lie).

And you should never (well, almost never) do x=(A'*A)\(A'*b) to solve a least squares problem. The only time you should do that is if you know your data is very well conditioned (not this case), and if A\b is slower than you'd like. Cholesky on the normal equations can be faster than QR on A ... but of course that's useless if A'*A is horribly conditioned.

(The condition number of A'*A is the square of the condition of A, so if cond(A) is 1e8 you could hope to get 8 digits of accuracy in your result. But cond(A'*A) is 1e16, so you get nothing in the result.)

"Tim Davis" <davis@cise.ufl.edu> wrote in message <haa339$1tk$1@fred.mathworks.com>...

> And you should never (well, almost never) do x=(A'*A)\(A'*b) to solve a least squares problem.

Right, but what you could do is

x = (A'*A + L'*L) \ (A'*b)

where L is appropriate Tikhonov matrix. In the standard form L could be sqrt(lambda)*speye(n).

It is anyway dangerous not using any kind of regularization while solving system involving a ill-conditioned matrix (it is *not* the case of OP's problem, since the problem he had is simply due to a bad scaling).

Bruno

How about:

invX = X\eye(size(X));

Dave

On 4 Okt, 02:45, "leo nidas" <bleonida...@yahoo.gr> wrote:
> Hi there,
>
> I want to solve or b the : (X'X)*b=X'*y. I have some issues with the inverse of X'*X.

Find a book on elementary linear algebra, read it,
contemplate it, and then use the ages-old standard way
to handle these things:

b = inv(X'X)X'y

Compute the SVD of X:

X = USV where S = diagonal, U'U = UU' =I, V'V = VV'= I

Insert in the expression to find

b = V'inv(S)U'y.

Since S is diagonal, S^-1 is easily computed.

Now, the matlab function PINV does all this,
so if you want a canned routine you can use it.
However, by writing this out in full, you can get
a useful solution by inspecting the elements on
the diagonal of S:

If |s(n,n)| ~ 0 for some n, the corresponding
term does not contribute in the data model.
However, the corresponding term 1/s(n,n) totally
dominates the pseudo inverse.

So in the pseudo inverse you will want to skip
all terms with singular values smaller than some
tolerance.

PINV does this, with the TOL argument that you
as user have to provide in advance - I suspect this
is an absolute tolerance. I would prefer a relative
tolerance, like

RTOL = eps('double')*max(diag(S));

in which case you will have to roll your own.

Rune

"Bruno Luong" <b.luong@fogale.findmycountry> wrote in message <haa48p$in7$1@fred.mathworks.com>...
> "Tim Davis" <davis@cise.ufl.edu> wrote in message <haa339$1tk$1@fred.mathworks.com>...
>
> > And you should never (well, almost never) do x=(A'*A)\(A'*b) to solve a least squares problem.
>
> Right, but what you could do is
>
> x = (A'*A + L'*L) \ (A'*b)
>
> where L is appropriate Tikhonov matrix. In the standard form L could be sqrt(lambda)*speye(n).
>
> It is anyway dangerous not using any kind of regularization while solving system involving a ill-conditioned matrix (it is *not* the case of OP's problem, since the problem he had is simply due to a bad scaling).
>
> Bruno

or redefine your initial system to include the L'L constraint and then to solve the redefined system
x=Aext\bext

Rune Allnor <allnor@tele.ntnu.no> wrote in message <6ed887e5-8dc1-4d0d-8d96-ef4ac3c7b5c7@j19g2000vbp.googlegroups.com>...
> On 4 Okt, 02:45, "leo nidas" <bleonida...@yahoo.gr> wrote:

>
> PINV does this, with the TOL argument that you
> as user have to provide in advance - I suspect this
> is an absolute tolerance. I would prefer a relative
> tolerance, like
>

Using PINV has two drawbacks:
1. It call for an expansive SVD
2. It does not work on sparse matrix

That's the reason I wrote a PSEUDOINVERSE on FEX based on QR factorization. Admittedly the factorization will be an expensive strategy for sparse matrix in term of storing space, and it is probably is impossible to apply on large scale problems (such as dimensions > 10e4, and yes there is a TONE of applications that need to solve such large system - forward or backward).

I might complete with an iterative regularized solver soon, even if it's no longer a "pseudo inverse" however both the practical scope of both have large overlapping.

> PINV does this, with the TOL argument that you
> as user have to provide in advance - I suspect this
> is an absolute tolerance.

Correct.

>I would prefer a relative
> tolerance, like
>
> RTOL = eps('double')*max(diag(S));

No big deal, just pass input TOL as the desired relative tolerance multiplied by norm/normest of the input matrix (which is the largest s(i,i)).

Bruno

"Petros Mpogiatzis" <painter@geo.auth.gr> wrote in message <haa8su$ige$1@fred.mathworks.com>...

>
> or redefine your initial system to include the L'L constraint and then to solve the redefined system
> x=Aext\bext

I don't see any obvious redefinition beside the obvious one:
Aext = A'*A + L'*L
bext = A'*b

Note that the Tikhonov regularization above has another equivalent formulation:

minimizing || A*x - b ||
such that ||L*x|| <= delta

delta is "some" number

Bruno

On 4 Okt, 16:05, "Bruno Luong" <b.lu...@fogale.findmycountry> wrote:
> Rune Allnor <all...@tele.ntnu.no> wrote in message <6ed887e5-8dc1-4d0d-8d96-ef4ac3c7b...@j19g2000vbp.googlegroups.com>...
> > On 4 Okt, 02:45, "leo nidas" <bleonida...@yahoo.gr> wrote:
>
> > PINV does this, with the TOL argument that you
> > as user have to provide in advance - I suspect this
> > is an absolute tolerance. I would prefer a relative
> > tolerance, like
>
> Using PINV has two drawbacks:
> 1. It call for an expansive SVD

First come up with something that actually works,
only then criticize speed.

> 2. It does not work on sparse matrix

Were there any sparse matrices in this thread?

> > PINV does this, with the TOL argument that you
> > as user have to provide in advance - I suspect this
> > is an absolute tolerance.
>
> Correct.
>
> >I would prefer a relative
> > tolerance, like
>
> > RTOL = eps('double')*max(diag(S));
>
> No big deal, just pass input TOL as the desired relative tolerance multiplied by norm/normest of the input matrix (which is the largest s(i,i)).

...which requires an equally expensive SVD to compute the
norm before one calls the SVD to compute the PINV...

Get a grip!

Rune

Rune Allnor <allnor@tele.ntnu.no> wrote in message <0d69cf8a-2563-4458-b112-401524512d0d@p36g2000vbn.googlegroups.com>...

>
> ...which requires an equally expensive SVD to compute the
> norm before one calls the SVD to compute the PINV...
>

WRONG again, the norm estimation use Krylov space technique (Arnoldi iteration). With very few iterations on Rayleigh ratio will converge toward the norm. There is a difference between finding FEW eigenvalue (here once) and all of them.

Your turn to get a grip and open your old algebra book Rune.

Bruno

Rune Allnor <allnor@tele.ntnu.no> wrote in message <0d69cf8a-2563-4458-b112-401524512d0d@p36g2000vbn.googlegroups.com>...

>
> First come up with something that actually works,
> only then criticize speed.

Yes I did.

>
> > 2. It does not work on sparse matrix
>
> Were there any sparse matrices in this thread?

Actually originally this thread is NOT even about PINV. It's about bad scaling. My comments on your post is direct toward the *general* use of PINV.

The OP only need to rescale his data, and a simple backslash will do. Look at my earlier post and John's post.
 
Bruno