diff --git a/project/testing.ipynb b/project/testing.ipynb index ac57874..c6e6b41 100644 --- a/project/testing.ipynb +++ b/project/testing.ipynb @@ -14,16 +14,23 @@ { "cell_type": "code", "execution_count": 2, - "id": "426735f0-859a-4f7d-83cf-f52f464a2737", + "id": "afd34800-7b9a-4006-9849-071fa197d962", "metadata": {}, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "x = [5.810770900369532e-16, 0.4999999999999994]\n" + ] + }, { "data": { "text/plain": [ "3-element Vector{Float64}:\n", - " 2.8498844747360086\n", - " -2.404529848146403\n", - " -0.3104749426792832" + " 0.9999999999999993\n", + " 1.9999999999999993\n", + " 2.999999999999999" ] }, "execution_count": 2, @@ -32,24 +39,26 @@ } ], "source": [ - "A = [1. 2; 3 4; 5 6]\n", - "A = qrhous!(A)\n", - "A * [1., 2, 3]" + "using LinearAlgebra: \\, qr\n", + "A = [1. 2.; 3. 4.; 5. 6.]\n", + "(Q, R) = qr(A)\n", + "x = A \\ [1., 2, 3]\n", + "@show x\n", + "A * x" ] }, { "cell_type": "code", "execution_count": 3, - "id": "afd34800-7b9a-4006-9849-071fa197d962", + "id": "8b3d71f8-2f4f-44f0-8c97-dd9d78706163", "metadata": {}, "outputs": [ { "data": { "text/plain": [ - "3-element Vector{Float64}:\n", - " 2.8498844747360064\n", - " -2.404529848146404\n", - " -0.31047494267928366" + "2-element Vector{Float64}:\n", + " -1.5012955407352218e-16\n", + " 0.5000000000000001" ] }, "execution_count": 3, @@ -58,10 +67,34 @@ } ], "source": [ - "using LinearAlgebra: qr!\n", - "A = [1. 2.; 3. 4.; 5. 6.]\n", - "(Q, R) = qr!(A)\n", - "Q * [1., 2, 3]" + "A = [1. 2; 3. 4.; 5. 6.]\n", + "A = qrhous!(A)\n", + "b = [1., 2, 3]\n", + "x = A \\ b" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "1883485c-b026-4a8a-8c34-b7c74ad79868", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "3-element Vector{Float64}:\n", + " 0.9999999999999991\n", + " 1.9999999999999996\n", + " 2.9999999999999996" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "A * x" ] } ], diff --git a/project/thinQR.jl b/project/thinQR.jl index 3df71ec..31be4dd 100644 --- a/project/thinQR.jl +++ b/project/thinQR.jl @@ -1,9 +1,9 @@ module thinQR -import Base: size, show, getproperty, getfield, propertynames, * -using LinearAlgebra: norm, I, triu, diagm +import Base: size, show, getproperty, getfield, propertynames, \, * +using LinearAlgebra: norm, I, triu, diagm, ldiv! -export QRthin, qrhous!, qyhous +export QRthin, qrhous!, qyhous, qyhoust mutable struct QRthin{T <: Real} A::AbstractVecOrMat{T} @@ -43,6 +43,15 @@ function qyhous(A::QRthin{T}, y::AbstractArray{T}) where T return z end +function qyhoust(A::QRthin{T}, y::AbstractArray{T}) where T + m, n = size(A.A) + z = deepcopy(y) + for j ∈ 1:n + z[j:m] = z[j:m] - A.A[j:m, j] * (A.A[j:m, j]' * z[j:m]) + end + return z +end + function calculateQ(A::QRthin{T}) where T if A.AQ != nothing return A.AQ @@ -61,7 +70,7 @@ function calculateR(A::QRthin{T}) where T return A.AR end m, n = size(A.A) - A.AR = triu(A.A[1:n, :], 1) + diagm(A.d) + A.AR = vcat(triu(A.A[1:n, :], 1) + diagm(A.d), zeros(m-n, n)) return A.AR end @@ -87,8 +96,18 @@ Base.propertynames(A::QRthin, private::Bool=false) = (:R, :Q, (private ? fieldna Base.size(A::QRthin) = size(getfield(A, :A)) -function *(A::QRthin{T}, b::AbstractVector{T}) where {T} - return qyhous(A, b) +function (\)(A::QRthin{T}, b::AbstractVector{T}) where T + n, m = size(A) + v = qyhoust(A, b) + x = zeros(m) + for j ∈ m:-1:1 + x[j] = (v[j] - x[j+1:m]' * A.A[j, j+1:m]) * A.d[j]^-1 + end + return x +end + +function (*)(A::QRthin{T}, x::AbstractVecOrMat{T}) where T + return qyhous(A, (A.R * x)) end end # module \ No newline at end of file