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
>  

Reply via email to