An implementation of DIOM #63

amontoison · 2018-09-27T21:04:06Z

Little brother of DQGMRES.

Currently, the factorization is H is done without pivoting.
Like DQGMRES, we can avoid the storage of H and I found a way to preconditioned it too.
This algorithm needs 2n + 1 vectors of mem coefficients.

coveralls · 2018-09-27T21:15:54Z

Coverage increased (+0.02%) to 99.506% when pulling dfa69c7 on Xx-fractale-xX:diom into c221384 on JuliaSmoothOptimizers:master.

amontoison · 2018-09-27T21:16:13Z

I didn't passed tests on julia 0.7:
ERROR: LoadError: LoadError: UndefVarError: A_mul_B! not defined

I don't understand why, It works on my laptop...
I tested it with include(Krylov.jl") and using Main.Krylov.

dpo · 2018-09-27T22:09:51Z

Possibly because this package hasn't yet been upgraded to 0.7 and so doesn't pull the correct version of LinearOperators?! You're probably using a recent enough version on your laptop. First we need to upgrade to 0.7 (in another PR).

amontoison · 2018-09-28T04:17:10Z

I found the problem, @abel renamed the function A_mul_B! by mul! with the version 5.0 of LinearOperators.
A_mul_B! is used to preallocated LinearOperators in variants.jl.
If we want that tests passed on julia 0.7, we need to phase out julia 0.6 and modify variants.jl.

codecov-io · 2018-09-28T04:40:02Z

Codecov Report

Merging #63 into master will increase coverage by 0.02%.
The diff coverage is 100%.

@@            Coverage Diff            @@
##           master     JuliaLang/julia#63      +/-   ##
=========================================
+ Coverage   99.48%   99.5%   +0.02%     
=========================================
  Files          19      20       +1     
  Lines        1542    1619      +77     
=========================================
+ Hits         1534    1611      +77     
  Misses          8       8

Impacted Files	Coverage Δ
src/variants.jl	`100% <ø> (ø)`	⬆️
src/Krylov.jl	`100% <ø> (ø)`	⬆️
src/diom.jl	`100% <100%> (ø)`

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update c221384...dfa69c7. Read the comment docs.

abelsiqueira · 2018-09-28T12:34:43Z

A_mul_B and similar were removed from 1.0 link and mul! with dispatch is the new way to handle it. LinearOperators v0.4.1 still uses A_mul_B and will only work on 0.6/0.7.
There should be a last release limiting LinearOperators on REQUIRE, then dropping 0.6.

amontoison · 2018-09-28T13:21:49Z

Thanks @abel !
I limited LinearOperators 0.4.1 in the REQUIRE with a previous commit but during the continous integration it's still LinearOperators 0.5.0 that is installed on Julia 0.7.

amontoison · 2018-10-04T06:39:38Z

LU factorization with partial pivoting don't break if H is singular. Furthermore the process is more stable, the biggest pivot is chosen at each step.

I'm quite sure that this version was never implemented. In the article Practical use of some Krylov subspace methods for solving indefinite and nonsymmetric linear systems, Youcef Saad talked about this variant but he made a mistake (from my point of view) in the update of x. I finally found an other way to do it which seems correct for me.

dpo · 2018-10-19T13:09:50Z

test/test_diom.jl

@@ -0,0 +1,71 @@
+include("get_div_grad.jl")


DIOM and DQGMRES should really be tested on unsymmetric systems as well.

amontoison · 2018-10-19T15:13:02Z

I also updated the reference of the article Practical use of some krylov subspace methods for solving indefinite and nonsymmetric linear systems.

dpo · 2018-10-19T15:16:36Z

src/diom.jl

+
+"""Solve the consistent linear system Ax = b using direct incomplete orthogonalization method.
+
+DIOM is similar to SYMMLQ for nonsymmetric problems.


I don't think this is accurate. FOM is to the Arnoldi process as CG is to the Lanczos process. I think it would be more accurate to say that DIOM is similar to CG with partial reorthogonalization.

Ok, I will rewrite the docString.

dpo · 2018-10-19T15:17:04Z

src/diom.jl

+It's a more economical algorithm based upon the LU factorization with partial pivoting.
+
+In the particular case where A is symmetric indefinite, DIOM is theorecally equivalent 
+to SYMMLQ but slightly more economical. 


I don't think this is true either; DIOM and SYMMLQ do not minimize the same quantities.

Yes, it's not norm(x) that is minimized. I think we could delete this sentence?

dpo · 2018-10-19T15:18:09Z

src/diom.jl

+    next_pos = mod(iter, mem) + 1 # Position corresponding to vₘ₊₁ in the circular stack V.
+
+    # Incomplete Arnoldi procedure.
+    z = M * V[:,pos] # Forms pₘ


I wonder why Github isn't rendering UTF-8 characters correctly. Is your editor set to use UTF-8 encoding?

Yes, my editor is set to use UTF-8 encoding. Github is rendering UTF-8 correctly for me, maybe it's your web browser?

Both Safari and Chrome seem to be missing characters...

dpo · 2018-10-19T15:18:54Z

Great, thank you!

dpo · 2018-10-20T13:35:51Z

src/diom.jl

+  # m-i+2 represents the indice of the diagonal where hᵢ.ₘ is located.
+  # In addition of that, the last column of Uₘ is stored in H.
+  L = zeros(mem) # Last mem Pivots of Lₘ.
+  p = zeros(Bool, mem) # Last mem permutations.


Maybe p could be a BitArray?

Yes, I didn't know that a Bool was an Uint8 in Julia and not a single bit.

dpo · 2018-10-20T13:37:41Z

src/diom.jl

+    end
+    # Compute hₘ₊₁.ₘ and vₘ₊₁.
+    H[1] = @knrm2(n, w) # hₘ₊₁.ₘ = ‖vₘ₊₁‖
+    if H[1] ≉ 0 # hₘ₊₁.ₘ ≈ 0 ⇒ "lucky breakdown"


Isn't there a tolerance involved in ≉? In my experience, we want to test for != 0.

Yes, there is a tolerance involved in ≈, x ≈ y return true if norm(x-y) ≤ rtol * max(|x|,|y|). With rtol that depend of the type of x and y. It's sqrt(eps) with eps the precision machine of Float16, Float32, Float64....

But when we do x ≈ 0, It's mean that |x| ≤ rtol * |x| ↔ x = 0. I used this notation because I find it nicer than ==. Also I like this notation because it means that x is a zero for the computer in the working precision but it's possible that it's not a zero in another precision or in exact arithmetic.

dpo · 2018-10-23T02:36:00Z

Thank you!

amontoison force-pushed the diom branch from ebad93c to a3e3f19 Compare September 27, 2018 21:08

amontoison force-pushed the diom branch from d30e959 to 37d0fdf Compare September 28, 2018 04:47

abelsiqueira closed this Oct 3, 2018

abelsiqueira reopened this Oct 3, 2018

amontoison added 6 commits October 14, 2018 18:52

An implementation of DIOM

d7fa15a

Correction of a comment about LU factorization in DIOM

776082f

LU factorization with partial pivoting in DIOM

ddffd4e

Add DIOM in test_variants.jl

e2e4f8e

Use exact residual norm in DIOM

1e5569b

Use a copy! instead of copy in DIOM

3ab0c65

amontoison force-pushed the diom branch from b74a0db to 3ab0c65 Compare October 14, 2018 22:52

abelsiqueira mentioned this pull request Oct 17, 2018

First Release #68

Closed

3 tasks

dpo reviewed Oct 19, 2018

View reviewed changes

test/test_diom.jl

@@ -0,0 +1,71 @@

include("get_div_grad.jl")

Copy link

Member

dpo Oct 19, 2018

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

DIOM and DQGMRES should really be tested on unsymmetric systems as well.

Tests with nonsymmetric matrix

20ca5f2

dpo reviewed Oct 19, 2018

View reviewed changes

Update the docstring of DIOM

147e280

dpo reviewed Oct 20, 2018

View reviewed changes

Use a BitArray instead of Array{Bool}

dfa69c7

dpo merged commit 075b1c6 into JuliaSmoothOptimizers:master Oct 23, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

An implementation of DIOM #63

An implementation of DIOM #63

amontoison commented Sep 27, 2018

coveralls commented Sep 27, 2018 •

edited

Loading

amontoison commented Sep 27, 2018

dpo commented Sep 27, 2018

amontoison commented Sep 28, 2018 •

edited

Loading

codecov-io commented Sep 28, 2018 •

edited by codecov bot

Loading

abelsiqueira commented Sep 28, 2018

amontoison commented Sep 28, 2018

amontoison commented Oct 4, 2018

dpo Oct 19, 2018

amontoison commented Oct 19, 2018

dpo Oct 19, 2018

amontoison Oct 19, 2018

dpo Oct 19, 2018

amontoison Oct 19, 2018

dpo Oct 19, 2018

amontoison Oct 19, 2018

dpo Oct 19, 2018

dpo commented Oct 19, 2018

dpo Oct 20, 2018

amontoison Oct 22, 2018

dpo Oct 20, 2018

amontoison Oct 22, 2018 •

edited

Loading

dpo Oct 23, 2018

dpo commented Oct 23, 2018


		"""Solve the consistent linear system Ax = b using direct incomplete orthogonalization method.

		DIOM is similar to SYMMLQ for nonsymmetric problems.

An implementation of DIOM #63

An implementation of DIOM #63

Conversation

amontoison commented Sep 27, 2018

coveralls commented Sep 27, 2018 • edited Loading

amontoison commented Sep 27, 2018

dpo commented Sep 27, 2018

amontoison commented Sep 28, 2018 • edited Loading

codecov-io commented Sep 28, 2018 • edited by codecov bot Loading

Codecov Report

abelsiqueira commented Sep 28, 2018

amontoison commented Sep 28, 2018

amontoison commented Oct 4, 2018

Choose a reason for hiding this comment

amontoison commented Oct 19, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpo commented Oct 19, 2018

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

amontoison Oct 22, 2018 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dpo commented Oct 23, 2018

coveralls commented Sep 27, 2018 •

edited

Loading

amontoison commented Sep 28, 2018 •

edited

Loading

codecov-io commented Sep 28, 2018 •

edited by codecov bot

Loading

amontoison Oct 22, 2018 •

edited

Loading