Solution of overdetermined linear systems in Complex Bigfloat erroneous

[andreasvarga]

The following is an example of an overdetermined linear system Ax = B with the following BigFloat complex data

julia> A = Complex{BigFloat}[1+im 1-im]
julia> B = Complex{BigFloat}[3+im]

The computed solution is

julia> x = A\B
2-element Vector{Complex{BigFloat}}:
   1.000000000000000000000000000000000000000000000000000000000000000000000000000017 - 0.5000000000000000000000000000000000000000000000000000000000000000000000000000086im
 -0.5000000000000000000000000000000000000000000000000000000000000000000000000000086 - 0.9999999999999999999999999999999999999999999999999999999999999999999999999999914im

The system being compatible, the residual must be zero, but it is not

julia> A*x-B
1-element Vector{Complex{BigFloat}}:
 -2.999999999999999999999999999999999999999999999999999999999999999999999999999965 - 0.9999999999999999999999999999999999999999999999999999999999999999999999999999741im

The correct solution, computed with ComplexF64 data, has the second element with changed sign:

julia> A1 = Complex{Float64}[1+im 1-im]
julia> B1 = Complex{Float64}[3+im]

julia> y = A1\B1
2-element Vector{ComplexF64}:
                1.0 - 0.4999999999999999im
 0.4999999999999998 + 1.0im
julia> A1*y-B1
1-element Vector{ComplexF64}:
 0.0 + 4.440892098500626e-16im

[oscardssmith]

This also appears with Complex{Float16} so I belive the problem is our fallback method.

[dkarrasch]

Must be our generic qrfactPivotedUnblocked then?

[aravindh-krishnamoorthy]

Hey, it seems I could quickly isolate the problem to here:

(Just copy the code below to "wide.jl" and run julia> include("wide.jl"))

using LinearAlgebra

function _wide_qr_ldiv!(A::QRPivoted{Complex{T}, Matrix{Complex{T}}, Vector{Complex{T}}, Vector{Int64}}, B::Matrix{Complex{T}}) where T
    m, n = size(A)
    minmn = min(m,n)
    mB, nB = size(B)
    lmul!(adjoint(A.Q), view(B, 1:m, :))
    R = A.R # makes a copy, used as a buffer below
    @inbounds begin
        if n > m # minimum norm solution
            τ = zeros(T,m)
            for k = m:-1:1 # Trapezoid to triangular by elementary operation
                x = view(R, k, [k; m + 1:n])
                τk = LinearAlgebra.reflector!(x)
                τ[k] = conj(τk)
                for i = 1:k - 1
                    vRi = R[i,k]
                    for j = m + 1:n
                        vRi += R[i,j]*x[j - m + 1]'
                    end
                    vRi *= τk
                    R[i,k] -= vRi
                    for j = m + 1:n
                        R[i,j] -= vRi*x[j - m + 1]
                    end
                end
            end
        end
        ldiv!(UpperTriangular(view(R, :, 1:minmn)), view(B, 1:minmn, :))
        if n > m # Apply elementary transformation to solution
            B[m + 1:mB,1:nB] .= zero(T)
            for j = 1:nB
                for k = 1:m
                    vBj = B[k,j]
                    for i = m + 1:n
                        vBj += B[i,j]*R[k,i]'
                    end
                    vBj *= τ[k]
                    B[k,j] -= vBj
                    for i = m + 1:n
                        B[i,j] -= R[k,i]*vBj
                    end
                end
            end
        end
    end
    return B
end

A = Complex{BigFloat}[1+im 1-im]
B = Complex{BigFloat}[3+im 0]
AF = Complex{Float64}[1+im 1-im]
BF = Complex{Float64}[3+im 0]
display(_wide_qr_ldiv!(qr(A,ColumnNorm()),reshape(B,2,1)))
display(_wide_qr_ldiv!(qr(AF,ColumnNorm()),reshape(BF,2,1)))

AF = Complex{Float64}[1+im 1-im]
BF = Complex{Float64}[3+im 0]
display(AF\BF)

c1 = rand()
c2 = rand()
d1 = rand()
C = Complex{BigFloat}[c1 c2]
D = Complex{BigFloat}[d1 0]
CF = Complex{Float64}[c1 c2]
DF = Complex{Float64}[d1 0]

display(_wide_qr_ldiv!(qr(C,ColumnNorm()),reshape(D,2,1)))
display(_wide_qr_ldiv!(qr(CF,ColumnNorm()),reshape(DF,2,1)))
CF = Complex{Float64}[c1 c2]
DF = Complex{Float64}[d1 0]
display(CF\DF)

Edit: new program shows that the algorithm in the above function fails also for Float64. It works otherwise because Float64 uses the BLAS+LAPACK routines.

Even more interesting with the random input, the results sometimes get reversed. I'll take a look at it later :)

[aravindh-krishnamoorthy]

The code above seems to have several problems. We solve it one by one.

1. The permutation matrix A.P is not applied.

Fix:

     # wrong code (in the previous comment)
     return B

     # correct code
     return  B[A.p,:] = B[1:n,:]

2. The application of the elementary transforms to the solution has wrong complex conjugation compared to LAPACK's zlarzb, which can be seen here: https://netlib.org/lapack/explore-html/d5/ddd/zlarzb_8f_source.html

Fix:

       if n > m # Apply elementary transformation to solution
            B[m + 1:mB,1:nB] .= zero(T)
            for j = 1:nB
                for k = 1:m
                    #  vBj = B[k,j]
                    vBj = B[k,j]'
                    for i = m + 1:n
                        # vBj += B[i,j]*R[k,i]'
                        vBj += B[i,j]'*R[k,i]'
                    end
                    vBj *= τ[k]
                    # B[k,j] -= vBj
                    B[k,j] -= vBj'
                    for i = m + 1:n
                        # B[i,j] -= R[k,i]*vBj
                        B[i,j] -= R[k,i]'*vBj'
                    end
                end
            end
        end

Commit is here: aravindh-krishnamoorthy/julia@53bf195

If someone experienced with this code reviews and agrees, I will submit a PR with the above commits (and some tests.)

Tests updated with 2x1, 3x1, and 3x2 dimensions for A, see attachment below. Download, rename to to "wide.jl", and run julia> include("wide.jl")

wide.jl.txt

[andreasvarga]

     # wrong code (in the previous comment)
     return B

     # correct code (e.g.)
     return transpose(A.P)*B

I think here it should be something similar to the previous code in qr.jl

     B[A.p,:] = B[1:n,:]
     return B

(not checked)

[aravindh-krishnamoorthy]

     # wrong code (in the previous comment)
     return B

     # correct code (e.g.)
     return transpose(A.P)*B

I think here it should be something similar to the previous code in qr.jl

     B[A.p,:] = B[1:n,:]
     return B

(not checked)

Thank you for this comment. Indeed transpose(A.P)*B was a typo in the comment. But, I compared the execution times for the "supposedly right" A.P*B and the proposal above based on your comment, i.e., B[A.p,:] = B[1:n,:]. As seen from below, the second method is better and included in the squashed commit: aravindh-krishnamoorthy/julia@53bf195. The comment above has also been updated.

The benchmark results are as follows (blank lines in the outputs have been stripped for presentation).

For A.P*B:

julia> include("wide.jl")
Original Test Case
2.482534153247272997077270316447523951720175849468581179858458680515081999974925e-16
0.0
2x1 Random Test Case
4.23814119446837277535136337493323174207783101841438353033438327862529691708443e-17
8.326672684688674e-17
3x1 Random Test Case
2.048541193521166664809196826486629262478121001976661594028921242391899921083791e-16
1.0007415106216802e-16
3x2 Random Test Case
2.792281291770672241662569322391667975812560099543899343836712833915874647116092e-16
1.5379475650849003e-16
julia> using BenchmarkTools
julia> @benchmark _wide_qr_ldiv!(qr(C,ColumnNorm()),reshape(DI,3,2))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  32.549 μs …  1.575 ms  ┊ GC (min … max): 0.00% … 97.27%
 Time  (median):     34.540 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   38.867 μs ± 57.187 μs  ┊ GC (mean ± σ):  6.25% ±  4.19%

  ▃▆██▆▃▂▁▁      ▁▁          ▁▁▁                              ▂
  ████████████▆▇█████▇█▇▇▇▇██████▆▅▆▆▆▅▄▅▅▄▅▃▄▃▂▄▄▅▃▂▅▄▄▄▂▂▄▂ █
  32.5 μs      Histogram: log(frequency) by time      69.3 μs <

 Memory estimate: 91.30 KiB, allocs estimate: 1814.

For B[A.p,:] = B[1:n,:]:

julia> include("wide.jl")
Original Test Case
2.482534153247272997077270316447523951720175849468581179858458680515081999974925e-16
0.0
2x1 Random Test Case
1.422769210168370174033308517599098978154697564002930685468287212995342972809115e-16
4.041272810440265e-16
3x1 Random Test Case
6.334071914299666997367508559030028378436801806234955251916140471763847712397077e-16
7.200419022730857e-16
3x2 Random Test Case
6.570734743074016968144327821288503092970420393144991406425972001689379925662375e-16
7.323081270985772e-16
julia> using BenchmarkTools
julia> @benchmark _wide_qr_ldiv!(qr(C,ColumnNorm()),reshape(DI,3,2))
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
 Range (min … max):  28.414 μs …  2.090 ms  ┊ GC (min … max): 0.00% … 95.35%
 Time  (median):     29.864 μs              ┊ GC (median):    0.00%
 Time  (mean ± σ):   33.650 μs ± 54.467 μs  ┊ GC (mean ± σ):  5.82% ±  3.60%

  ▄▇█▇▅▃▂▂▁▁   ▁▂▂▁▁      ▁▂▁▁                                ▂
  ████████████████████▇▆▇███████▇▇▇▆▇▆▆▄▆▃▆▅▄▅▄▅▅▅▅▅▄▄▃▄▅▅▁▃▅ █
  28.4 μs      Histogram: log(frequency) by time      58.4 μs <

 Memory estimate: 65.91 KiB, allocs estimate: 1317.

Furthermore, B[A.p,:] = B[1:n,:] is also better for single-column B (results not shown).

Revised commit: aravindh-krishnamoorthy/julia@53bf195

[aravindh-krishnamoorthy]

PR: JuliaLang/julia#49570

[dkarrasch]

Fixed by JuliaLang/julia#49570.