On Wednesday, April 2, 2014 10:46:10 AM UTC-5, Andreas Noack Jensen wrote:
>
> Maybe pivoted QR is best solution here.
>
It might be worthwhile trying out some other examples first. I can't
reproduce the Matlab answer from qrfact(A), qrfact(A;pivot=true) or
svdfact(A). That matrix does have a singular value of multiplicity 3 so
the last 3 columns are only defined up to orthogonal transformation.
Here's another example
julia> srand(1234321)
julia> A = randn((10,3))
10x3 Array{Float64,2}:
-0.160767 -1.1589 -2.05007
-1.48275 0.228077 -0.536018
-0.384419 -1.26447 0.920166
-1.3541 -0.0305313 -2.08457
1.03842 -1.73483 2.85367
-1.16268 1.10785 0.298823
-0.0202548 -0.0966725 0.30227
-0.747563 3.3424 -0.823113
1.32394 1.70391 1.65939
-1.26479 -1.17272 0.893986
julia> A = hcat(A, sum(A,2))
10x4 Array{Float64,2}:
-0.160767 -1.1589 -2.05007 -3.36974
-1.48275 0.228077 -0.536018 -1.79069
-0.384419 -1.26447 0.920166 -0.728726
-1.3541 -0.0305313 -2.08457 -3.46921
1.03842 -1.73483 2.85367 2.15726
-1.16268 1.10785 0.298823 0.243989
-0.0202548 -0.0966725 0.30227 0.185343
-0.747563 3.3424 -0.823113 1.77173
1.32394 1.70391 1.65939 4.68723
-1.26479 -1.17272 0.893986 -1.54353
julia> svdfact(A)
SVD{Float64,Float64}(10x4 Array{Float64,2}:
-0.445203 -0.0352878 -0.381188 0.452273
-0.230914 0.114977 0.403526 -0.592758
-0.0696868 -0.285529 0.263285 -0.042519
-0.46266 0.190456 0.0526681 -0.0314172
0.320501 -0.552839 0.141094 0.068795
0.0305325 0.194927 0.487114 0.633171
0.0282625 -0.0405129 0.0617822 0.0379557
0.191618 0.665113 0.208002 0.0473809
0.599258 0.136042 -0.0668694 0.0420263
-0.171123 -0.248158 0.5577 0.173786
,[8.96075,5.27466,2.59348,7.29607e-16],4x4 Array{Float64,2}:
0.248925 0.212351 0.39648 0.857757
-0.211592 0.824165 -0.516474 0.0960992
-0.802038 0.160183 0.571015 -0.0708406
-0.5 -0.5 -0.5 0.5 )
julia> qq = qrfact(A; pivot=true);
julia> qq[:R]
4x4 Array{Float64,2}:
7.70503 2.17425 2.20126 3.32952
0.0 -4.2384 1.47642 2.76199
0.0 0.0 1.87679 -1.87679
0.0 0.0 0.0 1.31664e-15
julia> qq[:Q]
10x10 QRPackedQ{Float64}:
-0.437343 0.0490778 0.388686 -0.419191
-0.232406 -0.173033 -0.381341 -0.115012
-0.094578 0.24982 -0.290425 0.665078
-0.450252 -0.22377 -0.0173727 -0.14866
0.279981 0.552938 -0.21007 -0.487088
0.0316663 -0.245138 -0.463803 -0.161661
0.0240548 0.0351485 -0.0666562 -0.040157
0.229944 -0.670641 -0.140445 -0.0845015
0.608334 -0.0899486 0.0626806 -0.158882
-0.200328 0.173925 -0.575773 -0.221894
I don't use Matlab. Could someone who does check what the result of
orth(A) is?
>
> 2014-04-02 11:39 GMT-04:00 Tomas Lycken <[email protected] <javascript:>
> >:
>
>> This is what I get from Matlab R2013a
>>
>> >> A
>>
>> A =
>>
>> 1 1 0 0 0
>> 1 1 0 0 0
>> 1 1 0 0 0
>> 1 1 0 0 0
>> 1 1 0 0 0
>> 1 0 1 0 0
>> 1 0 1 0 0
>> 1 0 1 0 0
>> 1 0 1 0 0
>> 1 0 1 0 0
>> 1 0 0 1 0
>> 1 0 0 1 0
>> 1 0 0 1 0
>> 1 0 0 1 0
>> 1 0 0 1 0
>> 1 0 0 0 1
>> 1 0 0 0 1
>> 1 0 0 0 1
>> 1 0 0 0 1
>> 1 0 0 0 1
>>
>> >> orth(A)
>>
>> ans =
>>
>> -0.2236 0.1060 0.3725 -0.0000
>> -0.2236 0.1060 0.3725 -0.0000
>> -0.2236 0.1060 0.3725 -0.0000
>> -0.2236 0.1060 0.3725 -0.0000
>> -0.2236 0.1060 0.3725 -0.0000
>> -0.2236 -0.3495 -0.0348 -0.1633
>> -0.2236 -0.3495 -0.0348 -0.1633
>> -0.2236 -0.3495 -0.0348 -0.1633
>> -0.2236 -0.3495 -0.0348 -0.1633
>> -0.2236 -0.3495 -0.0348 -0.1633
>> -0.2236 -0.0143 -0.1302 0.3645
>> -0.2236 -0.0143 -0.1302 0.3645
>> -0.2236 -0.0143 -0.1302 0.3645
>> -0.2236 -0.0143 -0.1302 0.3645
>> -0.2236 -0.0143 -0.1302 0.3645
>> -0.2236 0.2577 -0.2076 -0.2012
>> -0.2236 0.2577 -0.2076 -0.2012
>> -0.2236 0.2577 -0.2076 -0.2012
>> -0.2236 0.2577 -0.2076 -0.2012
>> -0.2236 0.2577 -0.2076 -0.2012
>>
>> On Wednesday, April 2, 2014 5:22:34 PM UTC+2, Douglas Bates wrote:
>>>
>>> On Tuesday, April 1, 2014 10:44:03 PM UTC-5, Gustavo Lacerda wrote:
>>>>
>>>> yes, I think that would be useful.
>>>>
>>>
>>> Do you know what Matlab's orth function returns for a rank deficient
>>> matrix? Suppose you ask for an orthogonal basis of the column space of a
>>> matrix like
>>>
>>> julia> hcat(fill(1.,(20,)), eye(4)[iceil([1:20] ./ 5),:])
>>> 20x5 Array{Float64,2}:
>>> 1.0 1.0 0.0 0.0 0.0
>>> 1.0 1.0 0.0 0.0 0.0
>>> 1.0 1.0 0.0 0.0 0.0
>>> 1.0 1.0 0.0 0.0 0.0
>>> 1.0 1.0 0.0 0.0 0.0
>>> 1.0 0.0 1.0 0.0 0.0
>>> 1.0 0.0 1.0 0.0 0.0
>>> 1.0 0.0 1.0 0.0 0.0
>>> 1.0 0.0 1.0 0.0 0.0
>>> 1.0 0.0 1.0 0.0 0.0
>>> 1.0 0.0 0.0 1.0 0.0
>>> 1.0 0.0 0.0 1.0 0.0
>>> 1.0 0.0 0.0 1.0 0.0
>>> 1.0 0.0 0.0 1.0 0.0
>>> 1.0 0.0 0.0 1.0 0.0
>>> 1.0 0.0 0.0 0.0 1.0
>>> 1.0 0.0 0.0 0.0 1.0
>>> 1.0 0.0 0.0 0.0 1.0
>>> 1.0 0.0 0.0 0.0 1.0
>>> 1.0 0.0 0.0 0.0 1.0
>>>
>>>
>>> Do you get 4 or 5 columns?
>>>
>>> The simple calculation gives you 5 columns but there should only be 4
>>> columns.
>>>
>>
>
>
> --
> Med venlig hilsen
>
> Andreas Noack Jensen
>