GMRES package?

[Benjamin Silbaugh]

I was thinking of implementing a GMRES solver in Julia for sparse linear systems. (I have some experience with this in C++.) However, before I do, I thought I would see if someone else has already done this, or is working on this? I looked through the package list for GMRES or sparse matrix solvers, but I didn't find anything.

The beauty of GMRES (and other Krylov methods) is that it obviates the need to explicitly formulate a (sparse) matrix. All you need is a function that computes the action of the linear operator (matrix) on a solution estimate. The disadvantage, however, is that can sometimes require a lot of memory--especially if you don't have a good pre-conditioner.

[John Myles White]

To my knowledge, no one has written any Krylov solvers so far for Julia. I'd be excited to see them.

-- John

[Amuthan Arunkumar Ramabathiran]

Hello Benjamin,

There was some discussion here on sparse iterative solvers for Julia a few months back and the consensus then was that it would be preferable to develop Julia wrappers to state-of-the-art software like PetSc or Trilinos rather than doing all the work again. I guess the problem is not about developing a simple Krylov solver library, but developing _efficient_ and _scalable_ ones. Having said that, I would personally be quite excited to see an iterative solver library entirely written in Julia. I was planning to work on this around December, but I would be very happy to contribute if you are interested in a collaborative work. 

Cheers,

Amuthan.

[Viral Shah]

I would love to see a julia Krylov solver package too. Given the language features, perhaps we can even do some interesting things with the API.

-viral

[Ben Sadeghi]

+1 for GMRES ... I'd find it very handy.

[ivan kissiov]

>>> The disadvantage, however, is that can sometimes require a lot of memory--especially if you don't have a good pre-conditioner.

     Benjamin, could you please explain how the quality of the pre-conditioner increases memory usage?

[Steven G. Johnson]

On Thursday, August 22, 2013 11:36:33 AM UTC-4, ivan kissiov wrote:

>>> The disadvantage, however, is that can sometimes require a lot of memory--especially if you don't have a good pre-conditioner.

     Benjamin, could you please explain how the quality of the pre-conditioner increases memory usage?

On an NxN linear system, GMRES with p iterations has O(Np) memory usage and requires O(Np^2) time, so it is unsustainable for large enough p if N is huge.    The solution is "restarted" GMRES, or GMRES(k), which restarts GMRES (in various ways) once p==k.  The problem is that restarted GMRES is no longer guaranteed to converge.

A preconditioner (a cheap approximate inverse of the matrix that is applied at each iteration) accelerates convergence of an iterative solver.   As a consequence, a good preconditioner allows you to use restarted GMRES(k) with a smaller k and still converge.   i.e. you can converge with lower memory usage as well as in less time.

--SGJ

[Benjamin Silbaugh]

Amuthan,

PetSc and Trilinos are impressive packages; however, my objectives (at least in the near term) were considerably less ambitious. I was thinking more along the lines of a basic/serial implementation of Krylov algorithms.

Of course, a parallel Krylov package would also be a nice contribution. Afterall, the future of scientific computing is clearly becoming more and more parallel. Whether or not parallelization at the level of a Krylov solver, which may be nested inside an "outer" solver/optimizer algorithm, is optimal for all applications and architectures, is another very interesting question.

[Hari Sundar]

The iterative methods in c++ library can be easily translated to Julia....

http://math.nist.gov/iml++/

I can start a package for this if enough people are interested..

Hari

[calu...@gmail.com]

I am also inclined to create a Krylov package entirely written in Julia. However after finding this post I would rather contribute than create a redundant package. What is the current status?

[Viral Shah]

There really isn't anything in place yet, as far as I know of. It would be great to start a KrylovSolvers.jl package.

-viral

[Benjamin Silbaugh]

Unfortunately, I haven't had much time to work on this. So, if you have the time to spare, go for it!

[calu...@gmail.com]

I will set up a GitHub repository this week and start coding by the end of next week.

Cheers!

[calu...@gmail.com]

I was a bit eager, so I have already coded the Conjugate Gradient algorithm. It haven't performed any optimisations yet. I will get to that after I have coded some more Krylov methods. Suggestions are welcome! Have a look at KrylovSolvers

[John Myles White]

Very, very cool. You might look into devectorizing your cg() code (manually or using the Devectorize.jl package) to make it a little faster and more memory conservative.

-- John

[calu...@gmail.com]

Thanks John for the tip. Memory usage is now reduced by 38% for a system matrix of size 2000x2000. Execution speed remained the same.

[Viral Shah]

I think that our language features should allow for a much nicer interface to the various solvers than matlab. The hard work is going to be the preconditioners and we will need to use stuff like petsc or other packages there.

-viral

[Tim Holy]

On a related topic, I have ported over the BSD-licensed LSQR from here:
http://www.stanford.edu/group/SOL/software/lsqr.html
but I haven't tested it seriously. Of interest? If so, do we want a "matrix
solvers" package? Or lots of little ones for particular algorithms? (The
latter would decrease load times...) Or is some of this suitable for base?

--Tim

[Viral Shah]

I think these interfaces and code should mature outside base for now. Eventually the common solvers and preconditioners should be in base.

-viral

[calu...@gmail.com]

@Viral: What do you exactly mean with nice interface? What kind of functionality would you like to have? And what in particular don't you like about interface of the MATLAB solvers?

[Viral Shah]

I am thinking of the kind of stuff we did with our Factorization objects. I don't have a specific idea on what we can do, but I intuitively feel that iterative solvers can also have some new API level thinking like we did with the Factorization objects.

-viral

[Tim Holy]

With these iterative solvers, I don't think the factorization is a good model,
since you don't ever compute a factorization.

Viral, we already have a very nice advantage over Matlab: AbstractMatrix. For
large-scale solvers, Matlab usually provides two versions, one that accepts an
ordinary matrix A and another where you supply function handles that compute
A*x and A'*x in cases where A is a linear operator that you don't want to
express explicitly as a matrix. (An example that comes to mind is convolution,
where it is better to use Fourier transforms than to create a whopping sparse
matrix.)

But in Julia all one has to do is define A_mul_B and A_mul_Bt for your type.
The solver can be written just in terms of A*x and A'*x, because these get
translated into your calls automatically. So all the core linear algebra code
can be written as it would look for matrices, even though it can handle
general linear operators just fine.

In other words, I think we get some wonderful Julian-esque features without
having to put a great deal of thought into the API of these solvers. It's nice
when good things happen for free.

--Tim

[Viral Shah]

I agree that we do not need factorization for the iterative solvers, but we will need to use the framework for factorization objects for preconditioners such as incomplete LU and incomplete Cholesky.

AbstractMatrix is certainly perfect here, and we get all that for free in julia.

The only other thing I believe that happens in matlab is that based on varargout, it does some more work. We will not have that flexibility, but I believe many of the results that are the output of iterative solvers in matlab are not expensive to compute.

-viral

[Jutho]

This looks very nice. I was wondering, would it not make sense to make the package more general (e.g. just call it Krylov) and also include sparse eigensolvers and other Krylov based methods, such as for computing the matrix exponential acting on a vector (which I happened to have implement just a few days ago). I would be happy to contribute to this effort, although I do not have much experience (but I have a copy of the newest edition of Golub and van Loan).

It would then make sense to first discuss a uniform style to be used throughout this package. For example, I personally prefer to optional parameters such as maxiter, top, krylovdim to be named arguments, so that one can specify only those that should deviate from the standard value.

[calu...@gmail.com]

I am first considering only sparse linear systems and later on I might add algorithms for sparse non-linear systems. However especially in the initial stages I would like to keep the scope of the package quite concise in order to focus on the main objective: solving sparse systems of equations efficiently. With that being said, sparse eigensolvers would be a great addition. However adding a function merely because it is loosely based on the same theory but does something completely different without directly contributing to the rest does not make a lot of sense.

I agree that the optional arguments should be changed to named arguments. I will change that next week.

[Jason Pries]

I need a fairly flexible GMRES implementation for the finite element code I'm developing for my PhD thesis. I've already started working on it so I can contribute what I have so far (real implementation) to the package.

[Carlos Baptista]

Jason, thanks for your GMRES implementation. As I already said on the pull request I will certainly check it out next week.

Currently I am in the final stages of my master thesis, however I have some plans of my own to develop a Finite Elements package for partitioned Fluid-Structure Interaction (FSI) simulations. That is also the main reason why I started the KrylovSolvers package and another package called PolyMath for polynomial-based arithmetic, interpolation, shape functions etc. Check it out!

[Benjamin Silbaugh]

I think this is a good start.

At some point you may want to consider generalizing your implementation to work with any linear operator; that is, a function that computes the action of a hypothetical matrix A on a vector x without having to explicitly construct the matrix A. Not having to explicitly construct a matrix is one of the major features of Krylov methods. I expect that this should be a straightforward revision, given Julia's support for first class functions. If someone still wishes to explicitly form a matrix A, then all they have to do is pass in a closure (a lamba function) that evaluates A*x.

If you decide to add support for nonlinear equations, then the linear operator may be a function that computes the directional derivative of some function f about some fixed point x0 in the direction x. In this case, use of closures/callbacks instead of sparse matrices becomes extremely useful. (See Brown and Saad's paper on Newton-Krylov methods.)

On Tuesday, October 8, 2013 9:18:42 PM UTC-4, Carlos Baptista wrote:

[Alan Edelman]

The basic versions of gmres,pcg, qmr,cgs  and the rest are really easy to write.

e.g. (gmres w/o restarts from the original paper)

[Benjamin Silbaugh]

@Tim,

What functions/operations are expected from a subclass of AbstractMatrix? (You're more experienced with Julia than I am.) If the expectation (or implicit policy) is that all subclasses of AbstractMatrix support more operators than just matrix multiplication, then I worry that this could cause issues with other Julia packages or client code that has structured it's methods around a stricter usage of the type hierarchy.

Perhaps the simplest solution here is to just pass a closure/callback function to the solver. Closures are super easy to construct in the REPL or driver script if needed, and new Julia users wouldn't need to fully understand the type hierarchy in order to use this package. Just my two cents...

[Jason Riedy]

And Alan Edelman writes:
> The basic versions of gmres,pcg, qmr,cgs and the rest are
> really easy to write. e.g. (gmres w/o restarts from the
> original paper)

In my experience so far, there will be Julia-specific oddities
relating to the types. I wish I had kept track of the little
bits rather than just fixing them, but sparse operations have
some limitations at the moment... Not huge, and many are obvious
(in hindsight) to people who know how type systems work.

If anyone developing iterative method packages for both dense
*and sparse* right-hand-sides could take notes along the way,
that will prove very useful. I didn't for something similar (a
specific splitting iterative method), and I regret it. You'll
run into more issues in the surrounding driver and test code than
the iterative routine itself.

[BTW, what is the preferred method for sending patches that does
*not* involve centralizing on github? I have far too many
accounts. My little patches aren't enough to justify getting
yet another account.]
--
Jason

[Stefan Karpinski]

On Fri, Oct 11, 2013 at 5:01 PM, Jason Riedy <ja...@lovesgoodfood.com> wrote:

[BTW, what is the preferred method for sending patches that does *not* involve centralizing on github?  I have far too many accounts. My little patches aren't enough to justify getting yet another account.]

You could use https://gist.github.com/ which lets you submit code and diffs without needing an account. However, I have to say that open source doesn't really work very well if people shift the burden of additional work (in this case creating a GitHub account and making a pull request, which GitHub allows you without leaving the web page of the file in question) onto the maintainers of the project. That said, we're always happy to accept fixes and contributions and if the change is small, cut-and-paste from an email will probably work.

[Jason Riedy]

And Stefan Karpinski writes:
> However, I have to say that open source doesn't really work

> very well if people shift the burden of additional work [...]

Nice try. Others here know me better. Also, see git itself (and
its original target project).
--
Jason

[Stefan Karpinski]

No, I actually checked out your site before writing that and saw that you are clearly a prolific open source contributor, but I don't think that really changes anything. Officially endorsing people sending Julia patches as gists or emails strikes me as a bad idea that will lead to mistakes and other chaos. In any case, I'm very happy to have you checking out Julia and would be even more thrilled to have you contribute.

[Steven G. Johnson]

On Friday, October 11, 2013 3:42:04 PM UTC-4, Benjamin Silbaugh wrote:

What functions/operations are expected from a subclass of AbstractMatrix? (You're more experienced with Julia than I am.) If the expectation (or implicit policy) is that all subclasses of AbstractMatrix support more operators than just matrix multiplication, then I worry that this could cause issues with other Julia packages or client code that has structured it's methods around a stricter usage of the type hierarchy.

Iterative methods like GMRES shouldn't be restricted to AbstractMatrices at all, because AbstractMatrices generally imply that the matrix elements are directly accessible.  Better not to declare a type of the argument at all ... allow the user to pass in any type that supports multiplication by a vector (i.e. use "duck typing").

The temptation to "over-type" functions is a continual temptation in Julia libraries.  (This is a bug in Julia's "eigs" function too, which is currently restricted to AbstractMatrix subtypes but shouldn't be.)

[Jutho]

Why should an AbstractMatrix imply that the matrix elements are available. You might be thinking of traditional sparse matrices, but a user is free to define any subtype of AbstractMatrix which his as fields all the variables that are required to reconstruct the matrix vector multiplication. I have added the following example to the testsuite of arpack recently:

import Base: size, issym, ishermitian

type CPM{T<:Base.LinAlg.BlasFloat}<:AbstractMatrix{T} # completely positive map

kraus::Array{T,3} # kraus operator representation

end

size(Phi::CPM)=(size(Phi.kraus,1)^2,size(Phi.kraus,3)^2)

issym(Phi::CPM)=false

ishermitian(Phi::CPM)=false

function *{T<:Base.LinAlg.BlasFloat}(Phi::CPM{T},rho::Vector{T})

rho=reshape(rho,(size(Phi.kraus,3),size(Phi.kraus,3)))

rho2=zeros(T,(size(Phi.kraus,1),size(Phi.kraus,1)))

for s=1:size(Phi.kraus,2)

As=slice(Phi.kraus,:,s,:)

rho2+=As*rho*As'

end

return reshape(rho2,(size(Phi.kraus,1)^2,))

end

The type CPM specifies a completely positive map, which is a kind of superoperator that maps positive definite matrices to positive definite matrices. It appears in the context of quantum information theory. Instead of directly specifying its matrix elements, I rather specify the Kraus operators:

http://en.wikipedia.org/wiki/Choi's_theorem_on_completely_positive_maps#Kraus_operators

And then I just define the matrix vector product using these Kraus operators. At no point do I directly compute the matrix elements of the CPM, since it would be much less efficient to construct the full matrix and then just multiply this one with the vector.

[Andreas Noack Jensen]

Unfortunately I have missed this thread and some double work has been done. I also needed some of this functionality in my thesis in which solve some large Toeplitz systems and I also settled on translating the MATLAB versions of the http://www.netlib.org/templates code instead of writing from scratch. I couldn't find any license information so I wrote an email to Jack Dongarra and it is all modified BSD. 

I think that we should try to work together on this and find the best interface. Mine is still pretty rough, but what I have done so far is here https://github.com/andreasnoackjensen/IterativeLinearSolvers.jl. I have tried to make the implementations slightly more memory efficient and therefore I have introduced the allocation free functions A_mul_B!(α,A,x,β,y) which tries to mimic BLAS' gemv and solve!(A,b) which writes the solution to b. The problem is that my solvers doesn't work with Julia as it is right now because of these new functions. There are some issues open on allocation free versions of multiplication and system solving but no decisions have been made as far as I know.

Regarding Steven G. Johnson's comments about AbstractMatrix, I think it would be useful to state more clearly how AbstractMatrix should be interpreted as I think there has been confusion about it more than once. 

[Viral B. Shah]

Yes, fully agree. Eigs actually does need to move away from AbstractMatrix.

-viral

[Viral B. Shah]

The templates code would have been my starting point as well, with the same optimizations you implemented.

There could be other types that support multiplication by a vector but are not AbstractMatrices. One could also just pass a function as in matlab. Requiring AbstractMatrix inputs just means that one has to write a bit more boilerplate code before using the iterative solvers or eigs with other types.

-viral

[Alan Edelman]

    I see students and others all the time whose world view is that matrices are collections of elements: dense or sparse; they haven't internalized the "linear operator" concept.

 

    So I show them how to solve (cumsum)(x)=b  and eigs(cumsum) .... and that rocks their world.

I point out that "cumsum" could be thought of tril(ones(n)) but the complexity would be silly.

In MATLAB, inv(cumsum)  is diff (but diff drops the 1st element):

>> A=@(x) cumsum(x);
>> b=(1:10)'.^2;
>> gmres(A,b)
gmres converged at iteration 9 to a solution with relative residual 8.8e-07.

ans =

    0.9999
    3.0002
    4.9999
    7.0001
    8.9999
   11.0000
   13.0000
   15.0000
   17.0000
   19.0000

>> diff(b)

ans =

     3
     5
     7
     9
    11
    13
    15
    17
    19

>> eigs(A,5)

ans =

     1
     1
     1
     1
     1

[Tim Holy]

Hah! Funny I didn't notice that my example of convolution which I used to
motivate this didn't provide elementwise access.

Basically I meant "abstraction". We really could use
https://github.com/JuliaLang/julia/issues/987
While I agree that overtyping is a problem, sometimes undertyping is too
because we rely on typing to determine which variant of a method gets called.

--Tim

[Alan Edelman]

In Julia backslash can take arguments, but I'm not sure if you can do

  A \("gmres") b

though I'm pretty sure you could do  ( \("gmres"))(A,b)

and  eig("arnoldi")(A)

[Jason Pries]

When I was working on my GMRES implementation, I noticed that calls to "dot" were 13x faster if I avoided sub-indexing by storing the subspace vectors as an "array of vectors" instead of in a matrix. Since most of the work in GMRES occurs constructing the orthonormal basis (for a large inner iterations) this speeds things up considerably. This will probably hold true for any method which constructs an orthonormal basis.

x    = rand(1000,2);
y    = Array(Array{Float64,1},2);
y[1] = rand(1000);
y[2] = rand(1000);
z    = rand(1000);

tic();
for i = 1:1e6
     a = dot(x,y[1]);
end
toc();

#elapsed time: 0.529180537 seconds

for i = 1:1e6
     a = dot(x,subs(z,1:1000,1));
end
toc();

#elapsed time: 6.843639198 seconds

In terms of interface, it makes sense to me to support supplying matrices and preconditioners as a functions. As a trivial example, I would expect both GMRES calls to be valid:

A    = rand(1000,1000);
L(x) = *(A,x);
b    = rand(1000);

x1 = gmres(A,b);
x2 = gmres(L,b);

The only real problem is with the different syntax (A*x, versus L(x)).

I strongly feel the interface should support both "left" and "right" preconditioners (where applicable). One thing that is a constant bother to me with the MATLAB interface (to GMRES at least) is that it only supports left preconditioners (i.e. solving   M*b = M*A*x), and not right preconditioners (i.e. b = A*M*y, x = M*y). The problem with left preconditioning is that instead of minimizing norm(A*x-b), norm(M*(A*x-b)) is minimized instead. A lot of the time, norm(M*(A*x-b)) being small does not necessarily translate to norm(A*x-b) being small. Right preconditioning avoids this problem.

I've tried to include these elements in my implementation https://github.com/priesj/KrylovSolvers.jl/blob/patch-1/src/solvers.jl but admittedly, I'm fairly new to Julia. I have a feeling I may be doing some superfluous things with types and/or breaking some conventions/best practices.

[John Myles White]

I believe the speed difference you're discovering is a result of current Julia (but not planned Julia) making a copy every time you construct a slice of an array.

-- John

[Stefan Karpinski]

Perhaps we need a LinearOperator abstraction above AbstractMatrix that doesn't imply efficient access to individual elements. Of course, given a linear operator, you can always compute `L[i,j] = e_i'*L*e_j`, but it won't necessarily be efficient, which changes the set of algorithms you may want to use.

[Jutho]

In the (physics) code I am working on, one of the first things I do is define AbstractLinearOperator as a typealias of AbstractMatrix. I can think of no example where linear solvers or eigensolvers should accept a kind of function that cannot be given the interpretation and cannot be implemented as a concrete subtype of AbstractMatrix. Any operation that is linear in its argument and maps a vector to a vector is by definition a linear operator, and for me AbstractMatrix was up till now Julia's name for that concept (which is why I am not sure whether we really need an additional layer above that, but I would not have anything against it either).

The only critique is indeed that Julia's type system does not have something like an interface or virtual methods to specify what the minimal set of methods is that a concrete subtype of AbstractMatrix should overload. Also, it would be better if the default answer for something like issym, isposdef or ishermitian of an AbstractMatrix is just false, rather then that these methods try to access all the individual matrix entries. So either the behavior of AbstractMatrix is changed so that it does not try to get the individual matrix elements, or an additional abstract type LinearOperator is introduced that supersedes AbstractMatrix and that has the aforementioned behavior.

But allowing linear solvers or eigensolvers to accept generic functions as argument, certainly sounds like a bad idea to me. You loose any possibility of encoding specific information of the AbstractMatrix/LinearOperator in the concrete type of the argument that you pass to the solver. This means that, like in matlab, you have to specify via additional arguments what the dimension of the vector space on which it acts, whether it is symmetric (or hermitian), etc. This does not fit within the philosophy of Julia's type system or of object-oriented design, where the solver in question should be able to extract these properties from the central argument (being the linear operator / abstract matrix) directly and the method dispatch system together with a well implemented user's type (which overloads the corresponding methods) makes everything work correctly. 

[Viral Shah]

Also, in order to get good performance, as we go ahead with iterative solvers, we need fast sparse matvec.

I am thinking of plugging in Rich Vuduc's OSKI library. Not sure if it is maintained or if there are replacements.

http://oski.sf.net

-viral

[Tim Holy]

On Saturday, October 12, 2013 10:26:58 AM Jason Pries wrote:
> When I was working on my GMRES implementation, I noticed that calls to
> "dot" were 13x faster if I avoided sub-indexing by storing the subspace
> vectors as an "array of vectors" instead of in a matrix. Since most of the
> work in GMRES occurs constructing the orthonormal basis (for a large inner
> iterations) this speeds things up considerably. This will probably hold
> true for any method which constructs an orthonormal basis.

In case it's helpful, there are some additional "hidden" versions of dot in
Base.LinAlg.BLAS. In particular, you can pass a pair of pointers and offsets,
which would (if it makes your life easier otherwise) allow you to use a matrix
to store the orthonormal basis.

To make it easier to use this, you could provide a wrapper

dot(A::AbstractMatrix, colA::Integer, B::AbstractMatrix, colB::Integer)

One could argue such a thing should be in base/; on the other hand, once we
have fast array views there will be no need for such a function.

--Tim

[Steven G. Johnson]

On Saturday, October 12, 2013 2:56:30 PM UTC-4, Stefan Karpinski wrote:

Perhaps we need a LinearOperator abstraction above AbstractMatrix

We can't put this "above" AbstractMatrix because AbstractMatrix is already a subtype of AbstractArray, and Julia doesn't support multiple inheritance.

My basic problem is that AbstractArray, if it means anything at all, means an indexable container.  i.e. an abstract matrix M should support M[i,j] if nothing else; while this can certainly be supported for arbitrary (finite-dimensional) linear operators (in a given basis), it may be extremely inefficient if the operator is stored implicitly.  

I really don't understand the need to declare a type here at all.  The *only* reason to declare a type in Julia is if you need to dispatch on it to distinguish it from other types.  Since eigs and GMRES are inherently only for linear operators, why must a matrix type be specified explicitly?  What's wrong with duck typing here?

--SGJ
 

[Steven G. Johnson]

On Saturday, October 12, 2013 7:27:11 PM UTC-4, Jutho wrote:

But allowing linear solvers or eigensolvers to accept generic functions as argument, certainly sounds like a bad idea to me. You loose any possibility of encoding specific information of the AbstractMatrix/LinearOperator in the concrete type of the argument that you pass to the solver. This means that, like in matlab, you have to specify via additional arguments what the dimension of the vector space on which it acts, whether it is symmetric (or hermitian), etc. This does not fit within the philosophy of Julia's type system or of object-oriented design, where the solver in question should be able to extract these properties from the central argument (being the linear operator / abstract matrix) directly and the method dispatch system together with a well implemented user's type (which overloads the corresponding methods) makes everything work correctly. 

I didn't say generic functions.  I said any type that supports * (or \ for inverse iteration).  You don't "lose" any flexibility this way.

1) You don't need to pass an extra argument for the size.  Just require the type to support size(A).

2) You can still have specialized versions of the function that are more efficient for specialized types, e.g. SymmetricTridiagonal.    (Note that declaring as AbstractMatrix does not help at all with this anyway, since AbstractMatrix provides no way to tell whether the matrix is. e.g. positive-definite).

--SGJ

 

[Steven G. Johnson]

On Sunday, October 13, 2013 9:08:04 AM UTC-4, Steven G. Johnson wrote:

I really don't understand the need to declare a type here at all.  The *only* reason to declare a type in Julia is if you need to dispatch on it to distinguish it from other types.  Since eigs and GMRES are inherently only for linear operators, why must a matrix type be specified explicitly?  What's wrong with duck typing here? 

To give a concrete example, insisting that the argument be an AbstractMatrix means that one cannot pass Factorization objects to eigs for efficient inverse-iteration Arnoldi, even those these objects represent perfectly good linear operators that support \ very efficiently, since Factorization is not a subtype of AbstractMatrix.

[Jutho]

On Sunday, October 13, 2013 3:08:04 PM UTC+2, Steven G. Johnson wrote:

On Saturday, October 12, 2013 2:56:30 PM UTC-4, Stefan Karpinski wrote:

Perhaps we need a LinearOperator abstraction above AbstractMatrix

We can't put this "above" AbstractMatrix because AbstractMatrix is already a subtype of AbstractArray, and Julia doesn't support multiple inheritance.

My basic problem is that AbstractArray, if it means anything at all, means an indexable container.  i.e. an abstract matrix M should support M[i,j] if nothing else; while this can certainly be supported for arbitrary (finite-dimensional) linear operators (in a given basis), it may be extremely inefficient if the operator is stored implicitly.  

These are indeed very good arguments (including those in your later posts such as Factorization objects). But maybe then, they seem to suggest, that a scientifically-aimed language such as Julia just needs two separated branches of array-like types in the type hierarchy. One for indexable containers, that support getindex, and that can have any type of elements, and has no specific relation to the mathemetical structure of a tensor but is just a useful programming construction, such as any array type in most other programming language.

And another one, which has the interpretation of a vector, linear operator, multilinear operator (tensor) depending on the dimension, where efficient access of the elements is not required, whose element type can probably be specified to be Numbers, and which have some other well-defined properties (such as the requirement that * between linear operator and vector is specified).

[Jutho]

The Factorization objects could then be subtypes of this second type of array.

It would just seem silly or confusing to explain to users: linear solvers / eigensolvers can be applied to any type that has the interpretation of a linear operator, that overloads the methods size, issym and ishermitian and implements the product with a vector,  but we don't have a common name/supertype for these kind of objects.

[Steven G. Johnson]

It's not at all silly to do duck typing rather than have a specific type. (Python uses duck typing extensively.) On the contrary, since Julia doesn't have multiple inheritance, it is essential to use duck typing where we can to maximize flexibility.

For example, Julia already uses "iterable" types extensively in its standard library, which is any type that implements next, start, and end functions rather than a type per se.

As another example, the quadgk integration function works on any type that supports +, -, * by scalar, and norm. That is, for any normed vector space...despite the fact that Julia does not have an abstract type for this.

[Jutho]

Ok, I guess you're right. I just had a hard time abandonding the idea of mimicking the hierarchy of mathematical/physical structures with Julia's type hierarchy. But I can find myself in your suggestion that linear solvers/eigensolvers should accept any type that have overloaded * with a vector, so that in the routines the matrix vector products can still be denoted as A*x and not as A(x) (if A were allowed to be an arbitrary function as suggested by others).

[Andreas Noack Jensen]

The project seems to be continued here http://bebop.cs.berkeley.edu/poski/. I tried Sparse BLAS in MKL which works right away with Julia's SparseMatrixCSC and there was a significant speed gain. I think it is based on the NIST Sparse BLAS so maybe that could be an alternative (when you have a working government as the NIST sites are shut down.)

2013/10/13 Viral Shah <vi...@mayin.org>

--
Med venlig hilsen

Andreas Noack Jensen

[Viral Shah]

I think the sparse BLAS in MKL is probably the easiest way to go - since we already have MKL licenses from Intel for including MKL in julia distributions.

-viral

[Viral Shah]

Andreas, can you commit your sparse BLAS interface to MKL?

-viral

[Andreas Noack Jensen]

Yes. I'll look into it.

2013/10/14 Viral Shah <vi...@mayin.org>

[Tony]

This is code that I used gmres function in Matlab.

function y = afun(x)

    y = [0; x(1:n-1)] + [x(2:n); 0];

end

n = 21;

x=ones(n,1);

b = afun(x);

v=gmres(@afun,-b)

but I do not know how to use it in julia. How to pass the afun into gmres when the input A (resp. Pl, Pr) is a function returning A*x. I tried it as below. Please help me!

function afun(x)

y = [0; x[1:n-1]] + [x[2:n]; 0];

return y  

end

n=21

x=ones(n,1)

b=afun(x)

v=IterativeSolvers.gmres(afun(x),-b)

Output: ERROR: LoadError: MethodError: `*` has no method matching *(::Array{Float64,1}, ::Array{Float64,1}) 

v=IterativeSolvers.gmres(afun,-b)

Output: ERROR: LoadError: MethodError: `one` has no method matching one(::Type{Any})

Thanks!

Vào 00:26:58 UTC+7 Chủ Nhật, ngày 13 tháng 10 năm 2013, Jason Pries đã viết: