Re: [julia-users] Re: higher rank sparse matrices
Perhaps https://github.com/JuliaComputing/NDSparseData.jl? --Tim On Sunday, August 21, 2016 8:14:48 AM CDT Kristof Cools wrote: > Just wondering whether there have emerged any packages for this in the > meantime. I need a rank 3 sparse matrix to implement a retarded potential > integral equation solver. The structure will have non zero entries for all > values of the first two indices and a varying but fixed length set for the > third index: > > A[i,j,k] != 0 for k0(i,j) <= k <= k1(i,j)
[julia-users] Re: higher rank sparse matrices
Just wondering whether there have emerged any packages for this in the meantime. I need a rank 3 sparse matrix to implement a retarded potential integral equation solver. The structure will have non zero entries for all values of the first two indices and a varying but fixed length set for the third index: A[i,j,k] != 0 for k0(i,j) <= k <= k1(i,j)
[julia-users] Re: higher rank sparse matrices
On Tuesday, September 15, 2015 at 11:02:26 PM UTC-4, Jack Poulson wrote: > > I believe that Tony is suggesting manually applying the sparse operator > rather than explicitly constructing it and then applying it. This is a > common (and significant) performance optimization when a sparse operator is > only used once or twice, as the construction of the sparse matrix is often > *more* expensive than applying it once. > > The answer to your question (since you are applying the sparse matrix from > the right to a row vector) would be a loop of the form: forall nonzero > triplets (i,j,value), y[j] += x[i] * value. > > On Tuesday, September 15, 2015 at 8:52:46 AM UTC-7, Frank Kampas wrote: >> >> Could you post a link to the part of the documentation that describes how >> to do that? >> >> On Tuesday, September 15, 2015 at 3:53:11 AM UTC-4, Tony Kelman wrote: >>> >>> Instead of constructing a sparse matrix in the inner loop it would be >>> more efficient to write an in place stencil kernel function to perform the >>> equivalent operation. >> >> This is what I'm doing now. It's a little faster. for i = 1:n for j = i:n if i == j lhs[rctr,i] = -2 * start[i] lhs[rctr,n+i] = -2 * start[n+i] lhs[rctr,2n+1] = -2 * (radii[i] - start[2n+1]) rhs[rctr] = -start[i]^2 - start[n+i]^2 + start[2n+1]^2- radii[i]^2 else lhs[rctr,i] = 2 * (start[i] - start[j]) lhs[rctr,j] = 2 * (start[j] - start[i]) lhs[rctr,n+i] = 2 * (start[n+i]-start[n+j]) lhs[rctr,n+j] = 2 *(start[n+j] - start[n+i]) rhs[rctr] = (start[i]-start[j])^2 + (start[n+i]-start[n+j])^2 + (radii[i]+radii[j])^2 end rctr += 1 end end
[julia-users] Re: higher rank sparse matrices
I believe that Tony is suggesting manually applying the sparse operator rather than explicitly constructing it and then applying it. This is a common (and significant) performance optimization when a sparse operator is only used once or twice, as the construction of the sparse matrix is often *more* expensive than applying it once. The answer to your question (since you are applying the sparse matrix from the right to a row vector) would be a loop of the form: forall nonzero triplets (i,j,value), y[j] += x[i] * value. On Tuesday, September 15, 2015 at 8:52:46 AM UTC-7, Frank Kampas wrote: > > Could you post a link to the part of the documentation that describes how > to do that? > > On Tuesday, September 15, 2015 at 3:53:11 AM UTC-4, Tony Kelman wrote: >> >> Instead of constructing a sparse matrix in the inner loop it would be >> more efficient to write an in place stencil kernel function to perform the >> equivalent operation. > >
[julia-users] Re: higher rank sparse matrices
Instead of constructing a sparse matrix in the inner loop it would be more efficient to write an in place stencil kernel function to perform the equivalent operation.
[julia-users] Re: higher rank sparse matrices
Could you post a link to the part of the documentation that describes how to do that? On Tuesday, September 15, 2015 at 3:53:11 AM UTC-4, Tony Kelman wrote: > > Instead of constructing a sparse matrix in the inner loop it would be more > efficient to write an in place stencil kernel function to perform the > equivalent operation.
[julia-users] Re: higher rank sparse matrices
I'm doing sequential linear programming on quadratic constraints. Using matrices makes this much more straight-forward. Without 4th rank matrices, I have to generate a large number of 2nd rank matrices for every iteration. However, I gather from your answer that only 2nd rank sparse matrices can be created. On Saturday, September 12, 2015 at 9:16:41 PM UTC-4, Tony Kelman wrote: > > In JuMP you can do indexing over constraints and variables with any number > of indexes. You probably don't need to worry about explicitly forming > constraint matrices at all, since the flattened individual indexes of > optimization variables and constraints are somewhat arbitrary and will > mostly likely be rearranged by presolve in any good optimization solver. > Just express the variables and constraints of the problem you want to solve > and JuMP will handle the messy transformations. > > > On Saturday, September 12, 2015 at 3:15:34 PM UTC-7, Frank Kampas wrote: >> >> >> >> On Saturday, September 12, 2015 at 12:09:11 PM UTC-4, Frank Kampas wrote: >>> >>> Is it possible to create sparse matrices with a rank other than 2? >>> >> >> I've been using 4th rank sparse matrices in Mathematica for circle >> packing. The constraints >> can be expressed using 2nd rank matrices and representing all the pairs >> of circles gives me >> another two dimensions. I'd like to move the code to Julia to take >> advantage of the access >> to various linear programming software packages. >> >
[julia-users] Re: higher rank sparse matrices
There aren't built-in data structures defined in Julia's standard library right now for higher-dimensional sparse matrices, no. But you can certainly come up with your own data structure and use it however you like. Are there any dimensions in your problem along which every 2-dimensional slice has the same sparsity pattern? If there are certain index sets you'll want to be doing reductions over more often than others, then you can pick a data structure appropriately. Or for maximum generality you could use a coordinate representation. On Sunday, September 13, 2015 at 7:28:42 AM UTC-7, Frank Kampas wrote: > > I'm doing sequential linear programming on quadratic constraints. Using > matrices makes this much more straight-forward. Without 4th rank matrices, > I have to generate a large number of 2nd rank matrices for every iteration. > However, I gather from your answer that only 2nd rank sparse matrices can > be created. > > On Saturday, September 12, 2015 at 9:16:41 PM UTC-4, Tony Kelman wrote: >> >> In JuMP you can do indexing over constraints and variables with any >> number of indexes. You probably don't need to worry about explicitly >> forming constraint matrices at all, since the flattened individual indexes >> of optimization variables and constraints are somewhat arbitrary and will >> mostly likely be rearranged by presolve in any good optimization solver. >> Just express the variables and constraints of the problem you want to solve >> and JuMP will handle the messy transformations. >> >> >> On Saturday, September 12, 2015 at 3:15:34 PM UTC-7, Frank Kampas wrote: >>> >>> >>> >>> On Saturday, September 12, 2015 at 12:09:11 PM UTC-4, Frank Kampas wrote: Is it possible to create sparse matrices with a rank other than 2? >>> >>> I've been using 4th rank sparse matrices in Mathematica for circle >>> packing. The constraints >>> can be expressed using 2nd rank matrices and representing all the pairs >>> of circles gives me >>> another two dimensions. I'd like to move the code to Julia to take >>> advantage of the access >>> to various linear programming software packages. >>> >>
[julia-users] Re: higher rank sparse matrices
This is the code that uses sparse matrices: for i = 1:n for j = i:n if i == j sp = start'*sparse([i,n+i,2n+1],[i,n+i,2n+1],[1,1,-1]) lhs[rctr,:] = -2*sp -2*radii[i]*c rhs[rctr] = -(sp*start)[1] -radii[i]^2 else sp = start'*sparse([i,j,n+i,n+j,i,j,n+i,n+j],[i,j,n+i,n+j,j,i,n+j,n+i],[1,1,1,1,-1,-1,-1,-1],2n+1,2n+1) lhs[rctr,:] = 2*sp rhs[rctr] = (sp*start)[1] + ( radii[i]+radii[j] )^2 end rctr += 1 end end With 4th rank sparse matrices, I could define them outside the loop. On Saturday, September 12, 2015 at 12:09:11 PM UTC-4, Frank Kampas wrote: > > Is it possible to create sparse matrices with a rank other than 2? >
[julia-users] Re: higher rank sparse matrices
Doesn't 'diag([1,1,1])` have a rank of 3? On Saturday, September 12, 2015 at 6:09:11 PM UTC+2, Frank Kampas wrote: > > Is it possible to create sparse matrices with a rank other than 2? >
[julia-users] Re: higher rank sparse matrices
On Saturday, September 12, 2015 at 12:09:11 PM UTC-4, Frank Kampas wrote: > > Is it possible to create sparse matrices with a rank other than 2? > I've been using 4th rank sparse matrices in Mathematica for circle packing. The constraints can be expressed using 2nd rank matrices and representing all the pairs of circles gives me another two dimensions. I'd like to move the code to Julia to take advantage of the access to various linear programming software packages.
[julia-users] Re: higher rank sparse matrices
In JuMP you can do indexing over constraints and variables with any number of indexes. You probably don't need to worry about explicitly forming constraint matrices at all, since the flattened individual indexes of optimization variables and constraints are somewhat arbitrary and will mostly likely be rearranged by presolve in any good optimization solver. Just express the variables and constraints of the problem you want to solve and JuMP will handle the messy transformations. On Saturday, September 12, 2015 at 3:15:34 PM UTC-7, Frank Kampas wrote: > > > > On Saturday, September 12, 2015 at 12:09:11 PM UTC-4, Frank Kampas wrote: >> >> Is it possible to create sparse matrices with a rank other than 2? >> > > I've been using 4th rank sparse matrices in Mathematica for circle > packing. The constraints > can be expressed using 2nd rank matrices and representing all the pairs of > circles gives me > another two dimensions. I'd like to move the code to Julia to take > advantage of the access > to various linear programming software packages. >
