In my case: to construct the generators of a Coxeter group from the Cartan matrix, you get the i-th one
by subtracting the ith line of the Cartan matrix from the corresponding line of the identity matrix.
If you have a Cartan matrix C, why not just do C - I to get the matrix whose rows are the generators?
Sorry, you did not get it. If the Cartan matrix C is nxn you get n generators which are nxn matrices.
The i-th matrix is the identity eye(C) modified as follows: you subtract the i-th line of C from the i-th line of
eye(C)
I see, so you are really computing G = I - eᵢ*eᵢ'*C = I - eᵢ*C[i,:]' where eᵢ is the unit coordinate vector eᵢ = zeros(n); eᵢ[i]=1. Yes, G = one(C); G[i,:] = C[i,:] (correction: G[i,:] .-= C[i,:]) is perhaps the most compact way to construct this.
(Assuming you need to construct it at all; depending on what you are doing, you should also consider working with G implicitly in matrix-free fashion: you can multiply a vector by G much more efficiently without constructing the matrix, given only C and i. Multiplying a vector by the full matrix G requires O(n²) operations, whereas multiplying by it implicitly using the definition requires only O(n) operations.)
In the above G[i,:] = C[i,:] should be G[i,:] -= C[i,:]
I construct other elements of the Coxeter group W by multiplying the generators in various ways — I am not interested only on the action on the vector space.
1 Like
You can probably still improve performance by not constructing the generator matrices explicitly. Multiplying m of your generator matrices via dense-matrix operations requires O(mn³) operations, but multiplying them implicitly requires only O(mn²) operations. Not constructing the product at all, but just storing it implicitly, you can compute the action of the product on a vector in O(mn) operations.
This seems like a good illustration of the general principle that explicitly constructing identity matrices is often a sign of suboptimal code.
2 Likes
Yes but in my applications n is small and m is large. For instance, the Weyl group of type E8 has
dimension n=8 but cardinality m about 6.0e8
Then you probably want to store the product explicitly (if you are going to act it on more than n vectors), but you could still gain a factor of ~8 by going from O(mn³) to O(mn²) to construct the n×n product of m generators by not computing the individual generator matrices explicitly.
Anyway, thank you for explaining the ultimate possible optimizations. I am currently trying to port large libraries
and ease of programming is my main concern. When my port succeeds it will be time to optimize.
That’s understandable. Just realize that if you are porting from e.g. Matlab-style code that constructs zillions (you mentioned 10⁸, yikes) of little temporary matrices, performance is likely to be quite poor in Julia until you rewrite it in a more Julian style. (Factors of 10–100 speedup from rewriting with Julia performance in mind are not uncommon.)
Check out GitHub - JuliaArrays/FillArrays.jl: Julia package for lazily representing matrices filled with a single entry which includes Eye(n,m).
3 Likes
Ah, nice. I have
as a utility in RigidBodyDynamics, a special case of Fill (a much better name than ConstVector). I might switch to using FillArrays.
Edit: I really think FillArrays this should be part of stdlib.
2 Likes
I think in the interest of keeping Stdlib as simple as possible, this is not necessary as users can easily add FillArrays.jl, and one would still need to type using FillArrays even if it were in Stdlib.
I actually think Fill, Ones, etc. should be the return types of fill, ones, etc. Some reasons:
- there are many cases where
fill,ones, etc. are used to construct a matrix that doesn’t need to be modified anymore, in which caseFill,Ones, etc. will result in much better performance. -
Base.fillis currently not very generic. It always returns anArray. Same forones,zeros, etc. Accordingly, there arespzeros,spdiagm,spones(on 0.6) methods for sparse arrays. This seems suboptimal. Wouldn’t it be better to have
SparseArray(zeros(3, 3))
Array(zeros(3, 3))
if you need mutable versions of a specific type?
I was under the impression that fill and ones were going to be deprecated in favour of something like Matrix(Rep(0), n, m). Though it looks like that may not happen anymore
Deprecate `ones`? · Issue #24444 · JuliaLang/julia · GitHub
Feel free to use any code from FillArrays.jl if you feel like making the proposed change as a PR.
PS I’m of the opinion that having synonyms like zeros(n,m) = Zeros(n,m) (or for that matter linspace(a,b,n) = LinSpace(a,b,n)) is not necessary: just use the capitalized version.
1 Like
Wouldn’t it be worthwhile to factor this thread into multiple (it appears at least two)?
It’s not perfect because the thread gradually diverged, but I think this is slightly better. Split the CAS discussion to Porting a CAS to Julia