weijianzhang/PolarFact.jl

A Julia package for the matrix polar decomposition
Switch branches/tags
Nothing to show
Fetching latest commit…
Cannot retrieve the latest commit at this time.

PolarFact

A Julia package for the matrix polar decomposition.

Install

To install the release version, type

`julia> Pkg.add("PolarFact")`

Overview

Every `m-by-n` matrix `A` has a polar decomposition `A=UH`, where the `m-by-n` matrix `U` has orthonormal columns if `m>n` or orthonormal rows if `m<n` and the `n-by-n` matrix `H` is Hermitian positive semidefinite. For a square matrix `A`, `H` is unique. If in addition, `A` is nonsingular, then `H` is positive definite and `U` is unique.

The polar decomposition is closely related to the singular value decomposition (SVD). In particular, if `A = P*S*Q'` is a singular value decomposition of A, then `U = P*Q'` and `H = Q*S*Q'` are the corresponding polar factors. The orthonormal polar factor `U` is the nearest orthonormal matrix to `A` in the Frobenius norm [1, Sec. 8.1].

[1] Nicholas J. Higham, Functions of Matrices: Theory and Computation, SIAM, Philadelphia, PA, USA, 2008.

This package provides the following algorithms for computing matrix polar decomposition:

• (Scaled) Newton's method

Reference: [2] Nicholas J. Higham, Computing the Polar Decomposition ---with Applications, SIAM J. Sci. Statist. Comput. Vol. 7, Num 4 (1986) pp. 1160-1174.

• the Newton Schulz method

This method can only apply to matrix `A` such that `norm(A) < sqrt(3)`.

Reference: [3] Günther Schulz, Iterative Berechnung der reziproken Matrix, Z. Angew. Math. Mech.,13:57-59, (1933) pp. 114, 181.

• a hybrid Newton method

Start with (scaled) Newton's method and switch to Newton-Schulz method when convergence is guaranteed.

Reference: [4] Nicholas J. Higham and Robert S. Schreiber, Fast Polar Decomposition of an arbitrary matrix, SIAM, J. Sci. Statist. Comput. Vol. 11, No. 4 (1990) pp. 648-655

• Halley's method

Reference: [5] Y. Nakatsukasa, Z. Bai and F. Gygi, Optimizing Halley's iteration for computing the matrix polar decomposition, SIAM, J. Mat. Anal. Vol. 31, Num 5 (2010) pp. 2700-2720.

• the QR-based Dynamically weighted Halley (QDWH) method [5]

• the SVD method

The scaled Newton iteration is a well known and effective method for computing the polar decomposition. It converges quadratically and is backward stable under the assumption that the matrix inverses are computed in a mixed backward forward stable way [6]. The QDWH is a cubic-rate convergent method. It is backward stable under the assumption that column pivoting and either row pivoting or row sorting are used in the QR factorization [6]. Without scaling, both type of methods can be slow when the matrix is ill-conditioned.

On many modern computers, matrix multiplication can be performed very efficiently. The Newton Schulz method requires two matrix multiplication while the (scaled) Newton method requires one matrix inversion. Thus the hybrid Newton is more efficient if matrix multiplication is 2 times faster than the matrix inversion [4].

Comparing to the SVD approach, the iterative algorithms are much more efficient when the matrix is nearly unitary (as in applications, for example, where a time-dependent matrix drifts from orthogonality due to rounding errors or other errors).

[6] Yuji Nakatsukasa and Nicholas J. Higham, Backward stability of iterations for computing the polar decomposition, SIAM, J. Matrix Anal. Appl. Vol. 33, No. 2, pp. 460-479.

Interface

The package provides a high-level function `polarfact`:

`	polarfact(A; alg, maxiter, tol, verbose)`

The meaning of the arguments:

• `A` : the input matrix of type `Matrix{T}`, where `T` could be `Float16`, `Float32`, `Float64`, `BigFloat`.

• `alg`: a symbol that indicates the factorization algorithm (default = `:newton`).

This argument accepts the following values:

• `:newton`: scaled Newton's method
• `:qdwh`: the QR-based Dynamically weighted Halley (QDWH) method
• `:halley`: Halley's method
• `:schulz`: the Newton Schulz method
• `:hybrid`: a hybrid Newton method
• `:svd`: the SVD method
• `maxiter`: maximum number of iterations (default = `100`).

• `tol` : tolerance (default = `cbrt(eps(T))`).

• `verbose` : whether to show procedural information (default = `false`), where `Iter` is the number of iterations, `Rel. err.` is equal to `norm(preU - U, Inf) / norm(preU, Inf)` and `Obj.` is equal to `norm(I - U'*U, Inf)`.

Note: `maxiter`, `tol` and `verbose` are not used for the SVD method.

The output has type `PolarFact.Result`, which is defined as

``````	immutable Result{T}
U::Matrix{T}               # unitary factor
H::Matrix{T}               # Hermitian positive semidefinite factor
niters::Union(Int, Nothing)      # number of iterations or Nothing
converged::Union(Bool, Nothing)  # whether the algorithm converges or Nothing
end
``````

Note: `niters` and `converged` are of type `Nothing` for the SVD method.

Examples

```julia> using PolarFact

julia> A = rand(6,6);

julia> r = polarfact(A);

julia> r.U
6x6 Array{Float64,2}:
0.78067    -0.0694445   0.470076    -0.292781   -0.0423338   0.277934
-0.0984991   0.200495   -0.350106    -0.167351    0.394071    0.802638
0.235931    0.376468    0.00966701   0.0259734   0.78062    -0.438716
0.38059    -0.347928   -0.256805     0.79899     0.105905    0.136195
0.21592     0.816347   -0.136636     0.261595   -0.445576    0.0362852
-0.365675    0.160542    0.756137     0.422827    0.154252    0.257271

julia> r.niters
6

julia> using MatrixDepot  # a test matrix collection

julia> A = matrixdepot("randsvd", 20, 10^15);  # test a very ill conditioned random matrix

julia> r = polarfact(A, alg = :newton, verbose = true);
Iter.    Rel. err.        Obj.
1     1.548278e+07     4.804543e+14
2     9.999002e-01     5.902800e+06
3     9.925194e-01     7.247316e+02
4     9.317962e-01     9.261149e+00
5     7.226745e-01     3.409441e-01
6     1.392861e-01     2.550612e-03
7     1.272341e-03     2.760504e-07
8     1.378506e-07     1.062337e-14

julia> r = polarfact(A, alg = :qdwh, verbose = true);
Iter.    Rel. err.        Obj.
1     1.018823e+00     2.294023e+00
2     9.492166e-01     2.113416e+00
3     6.766440e-01     7.363896e-01
4     1.337401e-01     8.038009e-04
5     1.208881e-04     4.126278e-13
6     6.184049e-14     2.242130e-15

julia> r = polarfact(A, alg = :halley);

julia> r.niters
34

8x8 Array{Float32,2}:
1.0   1.0   1.0   1.0   1.0   1.0   1.0   1.0
1.0  -1.0   1.0  -1.0   1.0  -1.0   1.0  -1.0
1.0   1.0  -1.0  -1.0   1.0   1.0  -1.0  -1.0
1.0  -1.0  -1.0   1.0   1.0  -1.0  -1.0   1.0
1.0   1.0   1.0   1.0  -1.0  -1.0  -1.0  -1.0
1.0  -1.0   1.0  -1.0  -1.0   1.0  -1.0   1.0
1.0   1.0  -1.0  -1.0  -1.0  -1.0   1.0   1.0
1.0  -1.0  -1.0   1.0  -1.0   1.0   1.0  -1.0

julia> r = polarfact(A)
Result{Float32}(8x8 Array{Float32,2}:
0.353553   0.353553   0.353553   0.353553  …   0.353553   0.353553   0.353553
0.353553  -0.353553   0.353553  -0.353553     -0.353553   0.353553  -0.353553
0.353553   0.353553  -0.353553  -0.353553      0.353553  -0.353553  -0.353553
0.353553  -0.353553  -0.353553   0.353553     -0.353553  -0.353553   0.353553
0.353553   0.353553   0.353553   0.353553     -0.353553  -0.353553  -0.353553
0.353553  -0.353553   0.353553  -0.353553  …   0.353553  -0.353553   0.353553
0.353553   0.353553  -0.353553  -0.353553     -0.353553   0.353553   0.353553
0.353553  -0.353553  -0.353553   0.353553      0.353553   0.353553  -0.353553,8x8 Array{Float32,2}:
2.82843     -1.49012e-8  -5.96046e-8  …  -1.49012e-8  -1.19209e-7  -1.49012e-8
-1.49012e-8   2.82843      5.96046e-8      5.96046e-8  -4.47035e-8  -7.45058e-8
-5.96046e-8   5.96046e-8   2.82843         0.0          8.9407e-8   -8.9407e-8
-7.45058e-8   1.49012e-8   2.98023e-8      1.3411e-7    1.3411e-7    1.49012e-7
-2.98023e-8   4.47035e-8   0.0             1.49012e-8   2.98023e-8  -7.45058e-8
-1.49012e-8   5.96046e-8   0.0         …   2.82843      1.63913e-7   7.45058e-8
-1.19209e-7  -4.47035e-8   8.9407e-8       1.63913e-7   2.82843      7.45058e-8
-1.49012e-8  -7.45058e-8  -8.9407e-8       7.45058e-8   7.45058e-8   2.82843   ,2,true)
```

Acknowledgements

The design of the package is inspired by NMF.jl.