I have just tried to build Julia with MKL to check this and with MKL I get the same result as in MATLAB. However, both results span the same space. With the MKL solution named U and the OpenBLAS solution named U2 I get
norm(U*U'-U2*U2') and the same is true for the pivoted QR. 2014-04-02 17:53 GMT-04:00 Douglas Bates <[email protected]>: > > 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]>: >> >>> 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 -- Med venlig hilsen Andreas Noack Jensen
