Merge pull request #47 from Marco-Congedo/dev

Dev

Project.toml

@@ -1,7 +1,7 @@
 name = "Diagonalizations"
 uuid = "9cd687f3-b62d-43f3-8fd3-ffcd9e581047"
 authors = ["Marco-Congedo <marco.congedo@gmail.com>"]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 CovarianceEstimation = "587fd27a-f159-11e8-2dae-1979310e6154"

README.md

@@ -1,6 +1,6 @@
 # Diagonalizations.jl
 
-| **Documentation**  | 
+| **Documentation**  |
 |:---------------------------------------:|
 | [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://Marco-Congedo.github.io/Diagonalizations.jl/dev) |
 
@@ -31,6 +31,10 @@ this package supports :
 - facilities to set the subspace dimension upon construction
 - diagonalization procedures for the case *m≥2* and *k≥2*.
 
+This package implements state-of-the-art **approximate joint diagonalization** algorithms. For some benchmarking see
+[here](https://github.com/Marco-Congedo/STUDIES/tree/master/AJD-Algos-Benchmark).
+
+
 ## Installation
 
 To install the package execute the following command in Julia's REPL:
@@ -40,7 +44,9 @@ To install the package execute the following command in Julia's REPL:
 ## Examples
 
 ```
-using Diagonalizations, Test
+
+using Diagonalizations, PosDefManifold, Test
+
 n, t=10, 100
 
 # generate an nxt data matrix
@@ -58,24 +64,25 @@ pC=pca(Symmetric((X*X')/t))
 
 # get only the first p eigenvectors, where p is the smallest integer
 # explaining at least 75% of the variance
-
-pX=pca(X; eVar=0.75) 
+pX=pca(X; eVar=0.75)
 
 Y=genDataMatrix(n, t)
 
 # maximum covariance analysis
-mXY=mca(X, Y) 
+mXY=mca(X, Y)
 
 # canonical correlation analysis
 cXY=cca(X, Y)
 
+# approximate joint diagonalization
+Xset=randP(5, 20)
+aXset=ajd(Xset; algorithm=:JADE)
+aXset=ajd(Xset; algorithm=:LogLike)
+
 # etc., etc.
 
 ```
 
-For the benchmarking of *Approxiate Joint Diagonalization* algorithms see 
-[here](https://github.com/Marco-Congedo/STUDIES/tree/master/AJD-Algos-Benchmark).
-
 ## About the Authors
 
 [Marco Congedo](https://sites.google.com/site/marcocongedo), is a research scientist of [CNRS](http://www.cnrs.fr/en) (Centre National de la Recherche Scientifique), working at [UGA](https://www.univ-grenoble-alpes.fr/english/) (University of Grenoble Alpes). **contact**: marco *dot* congedo *at* gmail *dot* com

docs/Project.toml

@@ -1,5 +1,5 @@
 authors = ["Marco-Congedo <marco.congedo@gmail.com>"]
-version = "0.1.7"
+version = "0.1.8"
 
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"

docs/src/Diagonalizations.md

@@ -265,7 +265,7 @@ The call
 p=pca(X; dims=2)
 ```
 
-uses the default SCM estimator and compute the PCA
+uses the default SCM estimator and computes the PCA
 from the ``100⋅100`` (rank-deficient) covariance matrix
 ``\frac{1}{30}XX^{T}``. The 'filter' `p.F` is in this case ``100⋅p``.
 
@@ -284,8 +284,8 @@ or `Hermitian` covariance matrices, the
 [mean](https://marco-congedo.github.io/PosDefManifold.jl/dev/riemannianGeometry/#Statistics.mean) function of package
 [PosDefManifold.jl](https://github.com/Marco-Congedo/PosDefManifold.jl)
 is used, since those matrices may be positive definite
-by construction, hence a mean using a metric acting on the Riemannian
-manifold of positive definite matrices may be used.
+by construction, hence a mean using a metric acting on the *Riemannian
+manifold of positive definite matrices* may be used.
 
 The constructors using this feature employ the following optional
 keyword arguments for regulating the computation of the mean:
@@ -462,15 +462,22 @@ upper case letters in order to allow reading in the REPL.
 - EEG: Electroencephalography
 - ERP: Event-Related Potentials
 - FOBI: Fourth-Order Blind Identification (Cardoso, 1989)
+- GAJD: Gauss AJD algorithm (Frobenius cost function), unpublished from the author of this package.
+- GLogLike: Gauss AJD algorithm (log-likelihood cost function), still experimental.
 - gCCA: generalized CCA
 - gMCA: generalized MCA
-- JADE: Joint Diagonalization of Eigenmatrices (Cardoso & Souloumiac, 1993)
+- JADE: Joint Diagonalization of Eigenmatrices AJD algorithm (Cardoso & Souloumiac, 1993)
+- JADEmax: JADE using the Riemannian gradient to order the rotations (Usevich et al., 2020)
+- LogLike: AJD algorithm optimizing the Log-Likelihood cost function (Pham, 2000)
+- LogLikeR: another implementation of LogLike for real data only.
 - LShrLW: Linear Shrinkage of Ledoit and Wolf (2004)
 - NShrLW: Non-linear Shrinkage of Ledoit and Wolf (2018)
+- mAJD: multiple AJD (jointly on several data sets)
 - MCA: Maximum Covariance Analysis
-- NoJoB: Non-orthogonal Joint BSS (Congedo et al., 2012)
-- OJoB: Orthogonal Joint BSS (Congedo et al., 2012)
+- NoJoB: Non-orthogonal mAJD algorithm (Congedo et al., 2012)
+- OJoB: Orthogonal mAJD algorithm (Congedo et al., 2012)
 - PCA: Principal Component Analysis (Pearson, 1901)
+- QNLogLike: Quasi-Newton LogLike (Ablin et al., 2019)
 - SCM: Sample Covariance Matrix
 - SOBI: Second-Order Blind Identification (Belouchrani et al., 1997)
 

docs/src/algorithms.md

@@ -14,15 +14,17 @@ Solvers differ from one another in several fashions:
 - their initialization
 - whether they support multi-threading or not
 
-To date seven solvers are implemented:
+To date nine solvers are implemented:
 
 - **Orthogonal joint blind source separation** (OJOB: Congedo et al., 2012[🎓](@ref))
 - **Non-orthogonal joint blind source separation** (NoJOB: Congedo et al., 2012[🎓](@ref))
-- **Log-likelyhood** (LogLike: Pham, 2001[🎓](@ref))
-- **Log-likelyhood Real** (LogLikeR: Pham, 2001[🎓](@ref))
-- **Quasi-Newton Log-likelyhood** (QNLogLike: Ablin et **al.**, 2019[🎓](@ref))
+- **Log-likelihood** (LogLike: Pham, 2001[🎓](@ref))
+- **Log-likelihood Real** (LogLikeR: Pham, 2001[🎓](@ref))
+- **Quasi-Newton Log-likelihood** (QNLogLike: Ablin et **al.**, 2019[🎓](@ref))
 - **Joint Approximate Diagonalization of Eigenmatrices** (JADE: Cardoso and Souloumiac, 1996[🎓](@ref))
+- **Joint Approximate Diagonalization of Eigenmatrices max** (JADEmax: Usevich, Li and Comon, 2020[🎓](@ref))
 - **Gauss Approximate Joint Diagonalization** (GAJD, unpublished, from the author)
+- **Gauss Log-Likelihood** (GLogLike, unpublished, from the author, still experimental)
 
 
 Their main characteristics and domain of application are summarized in the following table:
@@ -35,7 +37,10 @@ Their main characteristics and domain of application are summarized in the follo
 | LogLikeR | non-singular| no        | yes| log-likelihood | no | AJD |
 | QNLogLike| non-singular| no        | yes| log-likelihood | no | AJD |
 | JADE     | Orthogonal/Unitary| yes | no | least-squares | no | AJD |
+| JADEmax  | Orthogonal/Unitary| yes | no | least-squares | no | AJD |
 | GAJD     | non-singular| no        | no | least-squares | no | AJD |
+| GLogLike | non-singular| no        | yes | log-likelihood | no | AJD |
+
 
 
 ## using solvers
@@ -61,21 +66,12 @@ the actual approximate joint diagonalizer will be given by `init*B`, where `B` i
 `tol` is the tolerance for convergence of the solver.
 By default it is set to the square root of `Base.eps` of the nearest real type of the data input. If the solver encounters difficulties in converging, try setting `tol` in between 1e-6 and 1e-3.
 
-Argument `maxiter` is the maximum number of iterations allowed to the solver. The default values are given in the following table:
-
-| Algorithm  | max iterations |
-|:-----------|:---------------|
-| OJoB       | 1000 for real data, 3000 for complex data |
-| NoJoB      | 1000 for real data, 3000 for complex data |
-| Log-Like   | 120 for real data, 369 for complex data |
-| Log-LikeR  | 120 (real data only) |
-| QNLogLike  | 1000 (real data only) |
-| JADE       | 120 for real data, 360 for complex data |
-| GAJD       | 1000 (real data only) |
-
+Argument `maxiter` is the maximum number of iterations allowed to the solver.
+For all algorithms, by default `maxiter` is set to 1000 for real data and to
+3000 for complex data.
 
 If the maximum number of iteration
-is attained, a warning will be printed in the REPL.
+is attained, a warning is printed in the REPL.
 In this case, try increasing `maxiter` and/or `tol`.
 
 ## Multi-threading
@@ -88,10 +84,14 @@ mode paralellising several comptations over the dimension of the
 input matrices ``n``, if `threaded` is true, ``2n>x`` and ``x>1``,
 where ``x`` is the number of threads Julia is instructed to use.
 For QNLogLike, the algorithms run in multi-threaded
-mode paralellising several comptations over the number of matrices ``k``, if `threaded` is true, ``2k>x`` and ``x>1``.
+mode paralellising several comptations over the number of matrices ``k``,
+if `threaded` is true, ``2k>x`` and ``x>1``.
 
 Before running these methods you may want to set the number of threades
-Julia is instructed to use to the number of logical CPUs of your machine. Besides being optionally multi-threaded, OJoB and NoJoB algorithms heavely use BLAS. Before using these methods you may want to
-set `BLAS.set_num_threads(Sys.CPU_THREADS)` to the number of logical CPUs of your machine. If you are `using` any of the
+Julia is instructed to use to the number of logical CPUs of your machine.
+Besides being optionally multi-threaded, OJoB and NoJoB algorithms heavely use BLAS.
+Before using these methods you may want to
+set `BLAS.set_num_threads(Sys.CPU_THREADS)` to the number of logical CPUs
+of your machine. If you are `using` any of the
 package written by the author, all this is done automatically. See
 [these notes](https://marco-congedo.github.io/PosDefManifold.jl/dev/MainModule/#Threads-1).

docs/src/index.md

@@ -18,7 +18,7 @@ and run the examples you will find there.
 
 ### Disclaimer
 
-This version is throughoutly tested for the case of real data input, not tested at all for the case of complex data input.
+This version is throughoutly tested for both the case of real and complex data input.
 Independent reviewers for both the code and the documentation are welcome.
 
 
@@ -101,7 +101,7 @@ this package supports :
 
 ## Code units
 
-**Diagonalizations.jl** includes ten code units (.jl files):
+**Diagonalizations.jl** includes twelve code units (.jl files):
 
 | Unit   | Description |
 |:----------|:----------|
@@ -112,9 +112,11 @@ this package supports :
 | csp.jl | Unit implementing the CSP and CSTP|
 | ajd.jl | Unit implementing the AJD and the mAJD |
 | tools.jl | Unit containing general tools and internal functions|
+| Gajd.jl | Unit implementing the GAJD and GLogLike iterative algorithms (in the 'optim' folder)|
+| Jade.jl | Unit implementing the JADE and JADEmax iterative algorithms (in the 'optim' folder)|
 | JoB.jl | Unit implementing the OJoB and NoJoB iterative algorithms (in the 'optim' folder)|
-| LogLike.jl | Unit implementing the Log-Likelyhood iterative algorithms (in the 'optim' folder)|
-| Jade.jl | Unit implementing the Givens rotation iterative algorithms (in the 'optim' folder)|
+| LogLike.jl | Unit implementing the Log-Likelihood iterative algorithms (in the 'optim' folder)|
+| QnLogLike.jl | Unit implementing the quasi-Newton Log-Likelihood iterative algorithm (in the 'optim' folder)|
 
 Furthermore, all examples given at the end of the documentation of the filter
 constructors are collected as .jl units in the 'examples' folder.
@@ -184,9 +186,12 @@ P. Tichavsky, A. Yeredor (2009), [Fast approximate joint diagonalization
 incorporating weight matrices](https://staff.utia.cas.cz/files/Soutez_09/Tichavsky/WEDGE09.pdf), IEEE Trans. Signal Process., 57(3), 878–891.
 https://staff.utia.cas.cz/files/Soutez_09/Tichavsky/WEDGE09.pdf
 
-L. Tong, V. Soon, Y. Huang, R.W. Liu (1991). A necessary and sufficient condition Waveform-Preserving Blind
+L. Tong, V. Soon, Y. Huang, R.W. Liu (1991), A necessary and sufficient condition Waveform-Preserving Blind
 Estimation of Multiple Independent Sources. IEEE Transactions on Signal Processing, 41(7), 2461-2470.
 
+K. Usevich, J. Li, P. Comon (2020), [Approximate matrix and tensor diagonalization by unitary transformations: convergence of Jacobi-type
+algorithms](https://hal.archives-ouvertes.fr/hal-01998900v3/document),
+preprint: hal-01998900v3f.
 
 ## Contents