WIP: Hilbert/Banach space structures #27401
Conversation
|
Note: The new recursive definition of julia> using LinearAlgebra, StaticArrays
julia> x = SVector(1, 2); A = @SArray [1 2; 3 4];
julia> xx = [x, x]
2-element Array{SArray{Tuple{2},Int64,1,2},1}:
[1, 2]
[1, 2]
julia> AA = reshape([A, zero(A), zero(A), A], (2,2))
2×2 Array{SArray{Tuple{2,2},Int64,2,4},2}:
[1 2; 3 4] [0 0; 0 0]
[0 0; 0 0] [1 2; 3 4]
julia> norm(xx, 1)
6.0
julia> norm(reinterpret(eltype(x), xx), 1)
6.0while the old definition gives julia> norm(xx, 1)
4.47213595499958
julia> norm(reinterpret(eltype(x), xx), 1)
6.0 |
|
I have added a recursive definition of In order to get the old non-recursive definition of |
|
Currently, setting |
|
Much thanks for taking a shot at these changes @ranocha, and welcome to JuliaLang/julia! Superlative first contribution :). Apologies in advance for a bit of review latency: I hope to eke out review time over the next few evenings. Best! |
|
Thank you, @Sacha0 :) I won't possibly have as much time for this in the next days as I would like to have, but I will implement also some versions of I've fixed the error from |
|
I don't think we should have
|
As I think I wrote somewhere in #25565, I don't agree with this definition. Although it seems more natural at first glance, it is significantly less general — it requires that |
|
I think it is acceptable to have |
I'm very sorry, I did not remember such a comment in #25565. Thus, I thought that Edit: Should there be some possibility to get |
When setting
Could someone please give me a hint what has to be done to set |
It's easy to write this reduction yourself in cases where it is really needed. I think the default should be something that is as generic as possible, and in particular which doesn't break on arrays of custom types that only define
I'm guessing that the problem has to do with |
I have tried this on sunday with |
|
If it's working with |
It compiles well, but if I run I haven't seen this error before. I will try to see what's happening. |
| @@ -174,7 +174,7 @@ function (-)(J::UniformScaling, A::AbstractMatrix) | |||
| end | |||
|
|
|||
| inv(J::UniformScaling) = UniformScaling(inv(J.λ)) | |||
| norm(J::UniformScaling, p::Real=2) = abs(J.λ) | |||
| opnorm(J::UniformScaling, p::Real=2) = opnorm(J.λ, p) | |||
There was a problem hiding this comment.
Why not abs(J.λ)? This should be correct for all p since J.λ is a scalar, no? (Do we even have an opnorm method for scalars?)
There was a problem hiding this comment.
I thought someone might want to define some unusual number system with a strange norm as abs that is suited to the specific use case but does not yield a Banach algebra. But this might be rather unusual, so I could replace opnorm with abs. There is method opnorm(x::number, p::Real = 2) in https://github.com/JuliaLang/julia/pull/27401/files#diff-e5541a621163d78812e05b4ec9c33ef4R559.
There was a problem hiding this comment.
It is hard for me to imagine such an algebra, and norm(x::Number) defaults to abs(x) anyway, but I guess it doesn't matter as long as the opnorm method is defined, since performance is not critical here.
stdlib/LinearAlgebra/src/generic.jl
Outdated
|
|
||
| For any iterable containers `x` and `y` (including arrays of any dimension) of numbers (or | ||
| any element type for which `dot` is defined), compute the Euclidean dot product (the sum of | ||
| `dot(x[i],y[i])`) as if they were vectors. | ||
| any element type for which `inner` is defined), compute the inner product (or dot product |
There was a problem hiding this comment.
"Euclidean inner product" is perhaps more descriptive. (As opposed to some weighted inner product.)
Missing a close paren after "or scalar product".
There was a problem hiding this comment.
For me and the professors who have taught me linear algebra, Euclidean inner product implies that the corresponding field are the real numbers. Especially, if a vector space over the complex numbers is considered, the inner product has been called unitary inner product. But if Euclidean inner product means just an inner product in (american) English, I can of course rename it.
There was a problem hiding this comment.
The close parenthesis is after "or scalar product, i.e. the sum of inner(x[i],y[i]))". If it might be better to have this part of the additional explanation outside of the parentheses, I will of course change it.
There was a problem hiding this comment.
The reason for the qualifier "Euclidean" or similar is to distinguish it from some other weighted inner product. But since we write that it is the sum of inner(x[i],y[i]) explicitly, I guess we don't need any other qualifiers.
Yes, I would suggest having the i.e. the sum of ... part go outside of the parentheses, as this is necessary to precisely define what inner product we are using (as opposed to just a synonym).
There was a problem hiding this comment.
The docstring should also mention: dot(x,y) and x ⋅ y (where ⋅ can be typed by tab-completing \cdot in the REPL) are synonyms for inner(x,y).
It's also the first time I hear this argument I think. If the elements |
Yes? I'm not sure what you mean by "makes sense". The operation is well defined and is distinct from It is also not unreasonable to me — if you are using a vector of vectors as a data structure, it's not crazy to treat each subvector as a "unit" for purposes of e.g. an Linf norm. Of course, in cases where |
|
On Travis, the following error occurs: Locally, there have been no errors for me. I don't think I might have changed something that could cause the error above. Could someone please give me a hint? |
|
I'm seeing a similar |
Fair enough, but indeed, when it is just a vector of vectors, it is not compatible with the interpretation of this behaving like a block matrix (or block vector in this case). So I guess this need to be documented properly. With this argument, you also implicitly answered a question I posed in #25565: Is the specifier From your reasoning above, I conclude that the second behaviour is just fine. |
|
Thank you very much for your comments. I have updated the docstrings accordingly. |
|
There is a |
|
Should we just get rid of |
| @@ -529,15 +529,15 @@ end | |||
| diff(a::SparseMatrixCSC; dims::Integer) = dims==1 ? sparse_diff1(a) : sparse_diff2(a) | |||
|
|
|||
| ## norm and rank | |||
| vecnorm(A::SparseMatrixCSC, p::Real=2) = vecnorm(view(A.nzval, 1:nnz(A)), p) | |||
| norm(A::SparseMatrixCSC, p::Real=2) = norm(view(A.nzval, 1:nnz(A)), p) | |||
There was a problem hiding this comment.
do we actually need a view here?
There was a problem hiding this comment.
I think so… we don't want to make a copy of A.nzval[1:nnz(A)].
There was a problem hiding this comment.
but why not use A.nzval directly?
There was a problem hiding this comment.
The A.nzval might be longer than nnz(A).
There was a problem hiding this comment.
I see, thanks. Can that happen without manually messing with A.nzval?
There was a problem hiding this comment.
I thought it could happen with .= assignment to sparse matrices, but I don't remember clearly.
There was a problem hiding this comment.
Yes, the number of entries not matching the storage length can come about through e.g. broadcast.
|
@StefanKarpinski, it's not a huge deal to me, but I think it would be a shame to lose the infix notation and the familiar name to people coming from the physical sciences. On the other hand, I can see the argument for only having a single name for this function. |
ranocha
force-pushed
the
norm-opnorm-dot-inner
branch
from
eab4ff7
to
4282c4f
Compare
|
Thanks, @Sacha0. I've rebased on master. |
|
Awesome job @ranocha! |
|
Thanks! |
|
The change to julia> X1 = randn(3,2);
julia> X2 = [X1[i,:]' for i in 1:size(X1, 1)];
julia> X3 = reshape(X2, length(X2), 1);
julia> X1'X1
2×2 Array{Float64,2}:
2.58747 0.733045
0.733045 0.420716
julia> X2'X2
3.0081843634412744
julia> (X3'X3)[1]
2×2 Array{Float64,2}:
2.58747 0.733045
0.733045 0.420716 |
|
What do you propose as a way forward, @andreasnoack? |
|
We can't do anything before have we have versioned standard libraries or Julia 2.0 so, for now, it's just for the record. However, once we can adjust things, I'd prefer that |
|
I still think Ignorantly extrapolating my own coding style :-), maybe |
|
I just ran into this in the context of JuliaStats/StatsBase.jl#534 and JuliaStats/StatsBase.jl#518. I'm not sure if I caught all of the above (and related issues), but I'm curious - might there be a non-recursive variant of |
|
There is an open PR for such a function, but it has been mostly rejected because of doubts regarding naming and usefulness. |
|
Thanks for the link! |
This is based on the discussion in #25565 and is work in progress. It is my first pull request to julia. Thus, please be so kind to help me with improving this PR. Of course, it can also be used as a starting point for more experienced people.
The discussion in #25565 seems to have reached a consensus about the naming of
norm,vecnorm, and a newopnorm. Basically,normshould replace the oldvecnorm, since it is used in most cases to estimate some error or magnitude. Thus, avoiding the cost of the SVD for matrices seems to be good. Moreover, most users are probably not aware ofvecnormand just stick tonormsince it is okay in some cases. The old behaviour ofnorm(A::AbstractMatrix, p), i.e. the computation of the operator/matrix norm, is encoded inopnorm. Jutho has mentioned that he started implementing some similar behaviour in #25093. However, I can not see any commits or changed files there, so this is started from scratch.Moreover, the behaviour of genericnormfunctions is changed. Before, it has been along the lines ofSince it seems to be more natural (at least for me), I have changed it toThus, given the same coefficients, callingnormon a vector, a matrix or a vector of vectors yields the same result, which seems to be consistent in some sense.In order to keep
normanddot(or a newinner) consistent, I would also removevecdotand makedotrecursive, analogously tonorm, to give a behaviour asnorm(x) = sqrt(dot(x,x)). The discussion in #25565 about the naming ofdotdoes not seem to have reached a complete consensus. If nobody objected strongly, I would implement the basis inner/dot/scalar product asinnerand addconst dot = inneras proposed by Stefan Karpinski.This allows very generic code. One application I have in mind are multi-dimensional PDEs. For example, I would like to discretise a function in a cuboid as right-hand side for a Poisson problem. A straightforward way would be to use coefficients in an
b = Array{Float64,3}. Then, I have to solve a linear systemA \ b, whereAis a linear operator acting as negative Laplacian. Ifdotandnormwork as described above, I could simply use IterativeSolvers.jl for this task, sinceAis symmetric and positive definite.