diff --git a/NEWS.md b/NEWS.md index a46c5cccdfcd4..1597361af736f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -105,6 +105,9 @@ Deprecated or removed argument specifying the floating point type to which they apply. The old behaviour and `[get/set/with]_bigfloat_rounding` functions are deprecated ([#5007]) + * cholpfact and qrpfact are deprecated in favor of keyword arguments in + cholfact and qrfact. + [#4042]: https://github.com/JuliaLang/julia/issues/4042 [#5164]: https://github.com/JuliaLang/julia/issues/5164 [#5214]: https://github.com/JuliaLang/julia/issues/5214 @@ -135,6 +138,7 @@ Deprecated or removed [#5277]: https://github.com/JuliaLang/julia/issues/5277 [#987]: https://github.com/JuliaLang/julia/issues/987 [#2345]: https://github.com/JuliaLang/julia/issues/2345 +[#5330]: https://github.com/JuliaLang/julia/issues/5330 Julia v0.2.0 Release Notes ========================== diff --git a/base/deprecated.jl b/base/deprecated.jl index e28303b334be9..22191aba3642d 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -168,6 +168,14 @@ export PipeString # copy!(dest::AbstractArray, doffs::Integer, src::Integer) # in abstractarray.jl @deprecate copy!(dest::AbstractArray, src, doffs::Integer) copy!(dest, doffs, src) +@deprecate qrpfact!(A) qrfact!(A, pivot=true) +@deprecate qrpfact(A) qrfact(A, pivot=true) +@deprecate qrp(A, thin) qr(A, thin=thin, pivot=true) +@deprecate qrp(A) qr(A, pivot=true) +@deprecate cholpfact!(A) cholfact!(A, :U, pivot=true) +@deprecate cholpfact!(A,tol=tol) cholfact!(A, :U, pivot=true, tol=tol) +@deprecate cholpfact!(A,uplo,tol=tol) cholfact!(A, uplo, pivot=true, tol=tol) +@deprecate cholpfact(A) cholfact(A, :U, pivot=true) deprecated_ls() = run(`ls -l`) deprecated_ls(args::Cmd) = run(`ls -l $args`) diff --git a/base/linalg.jl b/base/linalg.jl index 8968705571b7d..89cd6da2c60ed 100644 --- a/base/linalg.jl +++ b/base/linalg.jl @@ -111,6 +111,7 @@ export svdvals!, svdvals, symmetrize!, + symmetrize_conj!, trace, transpose, tril, diff --git a/base/linalg/dense.jl b/base/linalg/dense.jl index df57d407783b8..c2863a3e1c6dc 100644 --- a/base/linalg/dense.jl +++ b/base/linalg/dense.jl @@ -423,7 +423,7 @@ function factorize!{T}(A::Matrix{T}) end return lufact!(A) end - return qrpfact!(A) + return qrfact!(A,pivot=true) end factorize(A::AbstractMatrix) = factorize!(copy(A)) @@ -443,7 +443,7 @@ function (\){T<:BlasFloat}(A::StridedMatrix{T}, B::StridedVecOrMat{T}) istriu(A) && return \(Triangular(A, :U),B) return \(lufact(A),B) end - return qrpfact(A)\B + return qrfact(A,pivot=true)\B end ## Moore-Penrose inverse diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index 299cc13084835..4209a27c58351 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -17,83 +17,51 @@ A_ldiv_B!(F::Factorization, B::StridedVecOrMat) = A_ldiv_B!(F, float(B)) Ac_ldiv_B!(F::Factorization, B::StridedVecOrMat) = Ac_ldiv_B!(F, float(B)) At_ldiv_B!(F::Factorization, B::StridedVecOrMat) = At_ldiv_B!(F, float(B)) +########################## +# Cholesky Factorization # +########################## type Cholesky{T<:BlasFloat} <: Factorization{T} UL::Matrix{T} uplo::Char end +type CholeskyPivoted{T<:BlasFloat} <: Factorization{T} + UL::Matrix{T} + uplo::Char + piv::Vector{BlasInt} + rank::BlasInt + tol::Real + info::BlasInt +end -function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol) +function cholfact!{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=-1.0) uplochar = string(uplo)[1] - C, info = LAPACK.potrf!(uplochar, A) - @assertposdef Cholesky(C, uplochar) info + if pivot + A, piv, rank, info = LAPACK.pstrf!(uplochar, A, tol) + return CholeskyPivoted{T}(A, uplochar, piv, rank, tol, info) + else + C, info = LAPACK.potrf!(uplochar, A) + return @assertposdef Cholesky(C, uplochar) info + end end cholfact!(A::StridedMatrix, args...) = cholfact!(float(A), args...) -cholfact!{T<:BlasFloat}(A::StridedMatrix{T}) = cholfact!(A, :U) -cholfact{T<:BlasFloat}(A::StridedMatrix{T}, args...) = cholfact!(copy(A), args...) -cholfact(A::StridedMatrix, args...) = cholfact!(float(A), args...) +cholfact{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Symbol=:U; pivot=false, tol=-1.0) = cholfact!(copy(A), uplo, pivot=pivot, tol=tol) +cholfact(A::StridedMatrix, uplo::Symbol=:U; pivot=false, tol=-1.0) = cholfact!(float(A), uplo, pivot=pivot, tol=tol) cholfact(x::Number) = @assertposdef Cholesky(fill(sqrt(x), 1, 1), :U) !(imag(x) == 0 && real(x) > 0) -chol(A::Union(Number, AbstractMatrix), uplo::Symbol) = cholfact(A, uplo)[uplo] -chol(A::Union(Number, AbstractMatrix)) = cholfact(A, :U)[:U] +chol(A::Union(Number, AbstractMatrix), uplo::Symbol=:U) = cholfact(A, uplo)[uplo] size(C::Cholesky) = size(C.UL) size(C::Cholesky,d::Integer) = size(C.UL,d) function getindex(C::Cholesky, d::Symbol) - C.uplo == 'U' ? triu!(C.UL) : tril!(C.UL) - if d == :U || d == :L - return symbol(C.uplo) == d ? C.UL : C.UL' - elseif d == :UL - return Triangular(C.UL, C.uplo) - end + d == :U && return triu!(symbol(C.uplo) == d ? C.UL : C.UL') + d == :L && return tril!(symbol(C.uplo) == d ? C.UL : C.UL') + d == :UL && return Triangular(C.UL, C.uplo) throw(KeyError(d)) end - -A_ldiv_B!{T<:BlasFloat}(C::Cholesky{T}, B::StridedVecOrMat{T}) = LAPACK.potrs!(C.uplo, C.UL, B) - -function det{T}(C::Cholesky{T}) - dd = one(T) - for i in 1:size(C.UL,1) dd *= abs2(C.UL[i,i]) end - dd -end - -function logdet{T}(C::Cholesky{T}) - dd = zero(T) - for i in 1:size(C.UL,1) dd += log(C.UL[i,i]) end - dd + dd # instead of 2.0dd which can change the type -end - -inv(C::Cholesky)=symmetrize_conj!(LAPACK.potri!(C.uplo, copy(C.UL)), C.uplo) - -## Pivoted Cholesky -type CholeskyPivoted{T<:BlasFloat} <: Factorization{T} - UL::Matrix{T} - uplo::Char - piv::Vector{BlasInt} - rank::BlasInt - tol::Real - info::BlasInt -end -function CholeskyPivoted{T<:BlasFloat}(A::StridedMatrix{T}, uplo::Char, tol::Real) - A, piv, rank, info = LAPACK.pstrf!(uplo, A, tol) - CholeskyPivoted{T}((uplo == 'U' ? triu! : tril!)(A), uplo, piv, rank, tol, info) -end - -chkfullrank(C::CholeskyPivoted) = C.rank 0 ? (dims < 3 ? size(A.hh, 1) : 1) : throw(BoundsError()) -function full{T<:BlasFloat}(A::QRPivotedQ{T}, thin::Bool=true) +function full{T<:BlasFloat}(A::QRPivotedQ{T}; thin::Bool=true) m, n = size(A.hh) Ahhpad = thin ? copy(A.hh) : [A.hh zeros(T, m, max(0, m - n))] LAPACK.orgqr!(Ahhpad, A.tau) end -print_matrix(io::IO, A::QRPivotedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, false), rows, cols) +print_matrix(io::IO, A::QRPivotedQ, rows::Integer, cols::Integer) = print_matrix(io, full(A, thin=false), rows, cols) ## Multiplication by Q from the Pivoted QR decomposition function *{T<:BlasFloat}(A::QRPivotedQ{T}, B::StridedVecOrMat{T}) diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index a4d449e14c025..1457282706643 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -1674,7 +1674,7 @@ for (posv, potrf, potri, potrs, pstrf, elty, rtyp) in (Ptr{BlasChar}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$rtyp}, Ptr{$rtyp}, Ptr{BlasInt}), &uplo, &n, A, &max(1,stride(A,2)), piv, rank, &tol, work, info) - @lapackerror + @assertargsok A, piv, rank[1], info[1] #Stored in PivotedCholesky end end diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index bacca28ffcc15..bbbd3d7980ea5 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -82,11 +82,9 @@ end size(A::Triangular, args...) = size(A.UL, args...) convert(::Type{Matrix}, A::Triangular) = full(A) -full(A::Triangular) = A.UL +full(A::Triangular) = A.uplo == 'U' ? triu!(A.UL) : tril!(A.UL) -getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? A.UL[i,j] : ((A.uplo == 'U') == (i < j) ? getindex(A.UL, i, j) : zero(T)) - -print_matrix(io::IO, A::Triangular, rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols) +getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? (A.unitdiag == 'U' ? one(T) : A.UL[i,j]) : ((A.uplo == 'U') == (i < j) ? getindex(A.UL, i, j) : zero(T)) istril(A::Triangular) = A.uplo == 'L' || istriu(A.UL) istriu(A::Triangular) = A.uplo == 'U' || istril(A.UL) diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 39187b527f396..defc0789a5f7c 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -59,56 +59,32 @@ Linear algebra functions in Julia are largely implemented by calling functions f Compute Cholesky factorization of a symmetric positive-definite matrix ``A`` and return the matrix ``F``. If ``LU`` is ``:L`` (Lower), ``A = L*L'``. If ``LU`` is ``:U`` (Upper), ``A = R'*R``. -.. function:: cholfact(A, [LU]) -> Cholesky +.. function:: cholfact(A, [LU,][pivot=false,][tol=-1.0]) -> Cholesky - Compute the Cholesky factorization of a dense symmetric positive-definite matrix ``A`` and return a ``Cholesky`` object. ``LU`` may be ``:L`` for using the lower part or ``:U`` for the upper part. The default is to use ``:U``. The triangular matrix can be obtained from the factorization ``F`` with: ``F[:L]`` and ``F[:U]``. The following functions are available for ``Cholesky`` objects: ``size``, ``\``, ``inv``, ``det``. A ``LAPACK.PosDefException`` error is thrown in case the matrix is not positive definite. + Compute the Cholesky factorization of a dense symmetric positive (semi)definite matrix ``A`` and return either a ``Cholesky`` if ``pivot=false`` or ``CholeskyPivoted`` if ``pivot=true``. ``LU`` may be ``:L`` for using the lower part or ``:U`` for the upper part. The default is to use ``:U``. The triangular matrix can be obtained from the factorization ``F`` with: ``F[:L]`` and ``F[:U]``. The following functions are available for ``Cholesky`` objects: ``size``, ``\``, ``inv``, ``det``. For ``CholeskyPivoted`` there is also defined a ``rank``. If ``pivot=false`` a ``PosDefException`` exception is thrown in case the matrix is not positive definite. The argument ``tol`` determines the tolerance for determining the rank. For negative values, the tolerance is the machine precision. .. function:: cholfact(A, [ll]) -> CholmodFactor Compute the sparse Cholesky factorization of a sparse matrix ``A``. If ``A`` is Hermitian its Cholesky factor is determined. If ``A`` is not Hermitian the Cholesky factor of ``A*A'`` is determined. A fill-reducing permutation is used. Methods for ``size``, ``solve``, ``\``, ``findn_nzs``, ``diag``, ``det`` and ``logdet``. One of the solve methods includes an integer argument that can be used to solve systems involving parts of the factorization only. The optional boolean argument, ``ll`` determines whether the factorization returned is of the ``A[p,p] = L*L'`` form, where ``L`` is lower triangular or ``A[p,p] = scale(L,D)*L'`` form where ``L`` is unit lower triangular and ``D`` is a non-negative vector. The default is LDL. -.. function:: cholfact!(A, [LU]) -> Cholesky +.. function:: cholfact!(A, [LU,][pivot=false,][tol=-1.0]) -> Cholesky ``cholfact!`` is the same as :func:`cholfact`, but saves space by overwriting the input A, instead of creating a copy. -.. function:: cholpfact(A, [LU]) -> CholeskyPivoted +.. function:: qr(A, [pivot=false,][thin=true]) -> Q, R, [p] - Compute the pivoted Cholesky factorization of a symmetric positive semi-definite matrix ``A`` and return a ``CholeskyPivoted`` object. ``LU`` may be ``:L`` for using the lower part or ``:U`` for the upper part. The default is to use ``:U``. The triangular factors contained in the factorization ``F`` can be obtained with ``F[:L]`` and ``F[:U]``, whereas the permutation can be obtained with ``F[:P]`` or ``F[:p]``. - The following functions are available for ``CholeskyPivoted`` objects: ``size``, ``\``, ``inv``, ``det``. - A ``LAPACK.RankDeficientException`` error is thrown in case the matrix is rank deficient. + Compute the (pivoted) QR factorization of ``A`` such that either ``A = Q*R`` or ``A[:,p] = Q*R``. Also see ``qrfact``. The default is to compute a thin factorization. Note that `R` is not extended with zeros when the full `Q` is requested. -.. function:: cholpfact!(A, [LU]) -> CholeskyPivoted +.. function:: qrfact(A,[pivot=false]) - ``cholpfact!`` is the same as ``cholpfact``, but saves space by overwriting the input A, instead of creating a copy. + Computes the QR factorization of ``A`` and returns either a ``QR`` type if ``pivot=false`` or ``QRPivoted`` type if ``pivot=true``. From a ``QR`` or ``QRPivoted`` factorization ``F``, an orthogonal matrix ``F[:Q]`` and a triangular matrix ``F[:R]`` can be extracted. For ``QRPivoted`` it is also posiible to extract the permutation vector ``F[:p]`` or matrix ``F[:P]``. + The following functions are available for the ``QR`` objects: ``size``, ``\``. When ``A`` is rectangular ``\`` will return a least squares solution and if the soultion is not unique, the one with smallest norm is returned. + The orthogonal matrix ``Q=F[:Q]`` is a ``QRPackedQ`` type when ``F`` is a ``QR`` and a ``QRPivotedQ`` then ``F`` is a ``QRPivoted``. Both have the ``*`` operator overloaded to support efficient multiplication by ``Q`` and ``Q'``. Multiplication with respect to either thin or full ``Q`` is allowed, i.e. both ``F[:Q]*F[:R]`` and ``F[:Q]*A`` are supported. A ``Q`` matrix can be converted into a regular matrix with ``full`` which has a named argument ``thin``. -.. function:: qr(A, [thin]) -> Q, R - - Compute the QR factorization of ``A`` such that ``A = Q*R``. Also see ``qrfact``. The default is to compute a thin factorization. Note that `R` is not extended with zeros when the full `Q` is requested. - -.. function:: qrfact(A) - - Computes the QR factorization of ``A`` and returns a ``QR`` type, which is a ``Factorization`` ``F`` consisting of an orthogonal matrix ``F[:Q]`` and a triangular matrix ``F[:R]``. - The following functions are available for ``QR`` objects: ``size``, ``\``. - The orthogonal matrix ``Q=F[:Q]`` is a ``QRPackedQ`` type which has the ``*`` operator overloaded to support efficient multiplication by ``Q`` and ``Q'``. Multiplication with respect to either thin or full ``Q`` is allowed, i.e. both ``F[:Q]*F[:R]`` and ``F[:Q]*A`` are supported. A ``QRPackedQ`` matrix can be converted into a regular matrix with ``full``. - -.. function:: qrfact!(A) +.. function:: qrfact!(A,[pivot=false]) ``qrfact!`` is the same as :func:`qrfact`, but saves space by overwriting the input A, instead of creating a copy. -.. function:: qrp(A, [thin]) -> Q, R, p - - Computes the QR factorization of ``A`` with pivoting, such that ``A[:,p] = Q*R``, Also see ``qrpfact``. The default is to compute a thin factorization. - -.. function:: qrpfact(A) -> QRPivoted - - Computes the QR factorization of ``A`` with pivoting and returns a ``QRPivoted`` object, which is a ``Factorization`` ``F`` consisting of an orthogonal matrix ``F[:Q]``, a triangular matrix ``F[:R]``, and a permutation ``F[:p]`` (or its matrix representation ``F[:P]``). - The following functions are available for ``QRPivoted`` objects: ``size``, ``\``. - The orthogonal matrix ``Q=F[:Q]`` is a ``QRPivotedQ`` type which has the ``*`` operator overloaded to support efficient multiplication by ``Q`` and ``Q'``. Multiplication with respect to either the thin or full ``Q`` is allowed, i.e. both ``F[:Q]*F[:R]`` and ``F[:Q]*A`` are supperted. A ``QRPivotedQ`` matrix can be converted into a regular matrix with ``full``. - -.. function:: qrpfact!(A) -> QRPivoted - - ``qrpfact!`` is the same as :func:`qrpfact`, but saves space by overwriting the input A, instead of creating a copy. - .. function:: bkfact(A) -> BunchKaufman Compute the Bunch Kaufman factorization of a real symmetric or complex Hermitian matrix ``A`` and return a ``BunchKaufman`` object. The following functions are available for ``BunchKaufman`` objects: ``size``, ``\``, ``inv``, ``issym``, ``ishermitian``. diff --git a/test/linalg.jl b/test/linalg.jl index 20dab6b92d619..8553d6bf2a63b 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -51,7 +51,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) @test_approx_eq l*l' apd # pivoted Choleksy decomposition - cpapd = cholpfact(apd) + cpapd = cholfact(apd, pivot=true) @test rank(cpapd) == n @test all(diff(diag(real(cpapd.UL))).<=0.) # diagonal should be non-increasing @test_approx_eq b apd * (cpapd\b) @@ -88,16 +88,16 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) # QR decomposition qra = qrfact(a) q,r = qra[:Q], qra[:R] - @test_approx_eq q'*full(q, false) eye(elty, n) - @test_approx_eq q*full(q, false)' eye(elty, n) + @test_approx_eq q'*full(q, thin=false) eye(elty, n) + @test_approx_eq q*full(q, thin=false)' eye(elty, n) @test_approx_eq q*r a @test_approx_eq a*(qra\b) b # (Automatic) Fat pivoted QR decomposition qrpa = factorize(a[1:5,:]) q,r,p = qrpa[:Q], qrpa[:R], qrpa[:p] - @test_approx_eq q'*full(q, false) eye(elty, 5) - @test_approx_eq q*full(q, false)' eye(elty, 5) + @test_approx_eq q'*full(q, thin=false) eye(elty, 5) + @test_approx_eq q*full(q, thin=false)' eye(elty, 5) @test_approx_eq q*r a[1:5,p] @test_approx_eq q*r[:,invperm(p)] a[1:5,:] @test_approx_eq a[1:5,:]*(qrpa\b[1:5]) b[1:5] @@ -105,8 +105,8 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) # (Automatic) Thin pivoted QR decomposition qrpa = factorize(a[:,1:5]) q,r,p = qrpa[:Q], qrpa[:R], qrpa[:p] - @test_approx_eq q'*full(q, false) eye(elty, n) - @test_approx_eq q*full(q, false)' eye(elty, n) + @test_approx_eq q'*full(q, thin=false) eye(elty, n) + @test_approx_eq q*full(q, thin=false)' eye(elty, n) @test_approx_eq q*r a[:,p] @test_approx_eq q*r[:,invperm(p)] a[:,1:5] @@ -125,7 +125,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) @test_approx_eq asym[1:5,1:5]*f[:vectors] scale(a610'a610*f[:vectors], f[:values]) @test_approx_eq f[:values] eigvals(asym[1:5,1:5], a610'a610) @test_approx_eq prod(f[:values]) prod(eigvals(asym[1:5,1:5]/(a610'a610))) - + # Non-symmetric generalized eigenproblem f = eigfact(a[1:5,1:5], a[6:10,6:10]) @test_approx_eq a[1:5,1:5]*f[:vectors] scale(a[6:10,6:10]*f[:vectors], f[:values]) @@ -149,7 +149,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) # singular value decomposition usv = svdfact(a) @test_approx_eq usv[:U]*scale(usv[:S],usv[:Vt]) a - + # Generalized svd gsvd = svdfact(a,a[1:5,:]) @test_approx_eq gsvd[:U]*gsvd[:D1]*gsvd[:R]*gsvd[:Q]' a @@ -163,7 +163,7 @@ for elty in (Float32, Float64, Complex64, Complex128, Int) x = tril(a)\b @test_approx_eq tril(a)*x b - + # Test null a15null = null(a[:,1:5]') @test rank([a[:,1:5] a15null]) == 10