Skip to content
This repository was archived by the owner on Nov 7, 2024. It is now read-only.

Jax gmres#707

Merged
chaserileyroberts merged 14 commits into
google:masterfrom
alewis:jax_gmres
Jul 8, 2020
Merged

Jax gmres#707
chaserileyroberts merged 14 commits into
google:masterfrom
alewis:jax_gmres

Conversation

@alewis

@alewis alewis commented Jul 2, 2020

Copy link
Copy Markdown
Contributor

Adds a Jax implementation of GMRES.

Comment thread tensornetwork/backends/jax/gmres.py Outdated
return x


@partial(jax.jit, static_argnums=(2,))

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Is this returning an arnoldi factorization?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

If so, do you think you can use the one that is implemented in the backend?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Well, the one in the backend returns the eigendecomposition; I need the V and H from the Krylov decomposition (actually I need V and the full QR decomposition of H). If Arnoldi were refactored to produce these it could be used.

However, a more efficient implementation, which I will add later, would also update the GMRES residual by a Givens rotation once per new Krylov vector. This becomes specialized enough that I think it makes more sense to keep things together with GMRES until proved otherwise, which is reflected in the code as written.

@mganahl mganahl Jul 3, 2020

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

sorry I was not precise in my wording. I meant the arnoldi factorization in jitted_functions.py. This one priduces exactly that. It' called _arnoldi_fact

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

would it make sense to move gmres into the jitted_functions.py file?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

My recent push should address the above

@chaserileyroberts chaserileyroberts left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

LGTM readability and design wise. Deferring to Martin's judgement for correctness.

@alewis

alewis commented Jul 7, 2020

Copy link
Copy Markdown
Contributor Author

The build works now

self._iram = jitted_functions._implicitly_restarted_arnoldi(libjax)
return self._iram(_CACHED_MATVECS[A], args, initial_state, num_krylov_vecs,
numeig, which, tol, maxiter)
imp_arnoldi = jitted_functions._implicitly_restarted_arnoldi(libjax)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

love that name!

"""
Helper function to generate jitted lanczos function used
in JaxBackend.eigsh_lanczos. The function `jax_lanczos`
Helper function to generate jitted lanczos function used

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Did you remove trailing white space on purpose, or is this from yapfing?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

I removed it on purpose with a vi script

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

why?

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

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

Trailing whitespace violates PEP8 and is in general Bad because for example it invisible modifies the behaviour of text editor commands to jump to the end of the line. Was there some reason to keep it?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

no, just curious. If you want it removed, then it would be great if you could prepare a PR that fixes this for the whole library (just run your script over all files). This way we'll have a cleaner code history

_, lambda x: norm_not_too_small)

return continue_iteration
initial_norm_typecaster = np.zeros((1,)) + eps + 1.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

you can do initial_norm = v.real.dtype.type(1.0+eps) instead

Given a linear mapping with (n x n) matrix representation
A = A_mv(*A_args) gmres_m solves
Ax = b (1)
where x and b are length-b vectors, using the method of

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

length-n vectors?


A_mv : A function `v0 = A_mv(v, *A_args, **A_kwargs)` where `v0` and
`v` have the same shape.
b : The `b` in `A @ x = b`.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

reverse order b and A_args

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

looks like A_args is missing from the docstring

return r, beta


#@partial(jax.jit, static_argnums=(2,))

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

why is the jit command commented?

Comment thread tensornetwork/backends/jax/jax_backend.py
backend.gmres(dummy_mv, b, A_kwargs=A_kwargs)


#jax_qr_dtypes = [np.float32, np.float64, np.complex64, np.complex128]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

remove commented code



#jax_qr_dtypes = [np.float32, np.float64, np.complex64, np.complex128]
jax_qr_dtypes = [np.float32]

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Also test for float64

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

and complex if possible

x, _ = backend.gmres(A_mv, b, x0=x0, num_krylov_vectors=n_kry, tol=tol)
solution = jax.numpy.array([2., 1.], dtype=dtype)
eps = jax.numpy.linalg.norm(jax.numpy.abs(solution) - jax.numpy.abs(x))
print(eps)

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

remove print

Comment thread tensornetwork/backends/jax/jax_backend_test.py
def A_mv(x):
return A @ x
b = A_mv(solution)
tol = b.size * jax.numpy.finfo(dtype).eps

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

good to know!

@alewis

alewis commented Jul 7, 2020

Copy link
Copy Markdown
Contributor Author

Martin's comments should be dealt with

@mganahl mganahl left a comment

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

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

Can you add the A_args to the docstring, then this is good to go

@mganahl

mganahl commented Jul 8, 2020

Copy link
Copy Markdown
Contributor

Hey @alewis, did you check that your modifications to eigsh_lanczos and eigs do not trigger unneccesary tracing? I checkout out your PR and on my laptop it seems it does

@alewis

alewis commented Jul 8, 2020

Copy link
Copy Markdown
Contributor Author

I find that tracing is triggered by the first run only, as in the screenshot. However, running the cell twice does trace twice; this also appears to have been true in the older code.

Screen Shot 2020-07-08 at 11 31 09 AM

@mganahl

mganahl commented Jul 8, 2020

Copy link
Copy Markdown
Contributor
import jax
import tensornetwork as tn
import numpy as np
import time
D=200
dtype=np.float32
matrix = jax.numpy.array(np.random.rand(D,D).astype(dtype))
vector = jax.numpy.array(np.random.rand(D,).astype(dtype))
@jax.jit
def matvec_jax_matrix(vec,matrix):
    return jax.numpy.tensordot(matrix, vec,([1],[0]))
jax_backend = tn.backends.jax.jax_backend.JaxBackend()
ncv=10
t1 = time.time()
eta_j, U_j = jax_backend.eigsh_lanczos(matvec_jax_matrix,[matrix],vector,num_krylov_vecs = ncv,numeig=1, 
                                       reorthogonalize=False)
print('jax eigvals:', eta_j)
t2 = time.time()
eta_j, U_j = jax_backend.eigsh_lanczos(matvec_jax_matrix,[matrix],vector,num_krylov_vecs = ncv,numeig=1, 
                                       reorthogonalize=False)

print('jax eigvals:', eta_j)
t3 = time.time()
print('jax first:', t2 - t1)
print('jax second:', t3 - t2)

could you run this and send the result?

@alewis

alewis commented Jul 8, 2020 via email

Copy link
Copy Markdown
Contributor Author

@mganahl

mganahl commented Jul 8, 2020

Copy link
Copy Markdown
Contributor

sorry can't see it

@mganahl

mganahl commented Jul 8, 2020

Copy link
Copy Markdown
Contributor

looks good, I'll pull it in as soon as building's finished

@chaserileyroberts chaserileyroberts merged commit db81bc0 into google:master Jul 8, 2020
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.

Labels

Projects

None yet

Development

Successfully merging this pull request may close these issues.

5 participants