vexilligera commented on a change in pull request #17014: [NumPy] Add NumPy support for norm URL: https://github.com/apache/incubator-mxnet/pull/17014#discussion_r368790115
########## File path: src/operator/numpy/linalg/np_norm-inl.h ########## @@ -0,0 +1,830 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one + * or more contributor license agreements. See the NOTICE file + * distributed with this work for additional information + * regarding copyright ownership. The ASF licenses this file + * to you under the Apache License, Version 2.0 (the + * "License"); you may not use this file except in compliance + * with the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, + * software distributed under the License is distributed on an + * "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY + * KIND, either express or implied. See the License for the + * specific language governing permissions and limitations + * under the License. + */ + +/*! + * Copyright (c) 2019 by Contributors + * \file np_norm-inl.h + * \brief norm + */ +#ifndef MXNET_OPERATOR_NUMPY_LINALG_NP_NORM_INL_H_ +#define MXNET_OPERATOR_NUMPY_LINALG_NP_NORM_INL_H_ + +#include <mxnet/operator_util.h> +#include <vector> +#include <limits> +#include <cmath> +#include "../../tensor/la_op.h" +#include "../../tensor/la_op-inl.h" +#include "../../tensor/init_op.h" +#include "./broadcast_reduce_op_customized.h" +#include "./np_gesvd-inl.h" +#include "../np_matrix_op-inl.h" + +namespace mxnet { +namespace op { + +namespace mshadow_op { +/*! \brief Lp-norm power reducer */ + +struct nrmlp { + double lp; + MSHADOW_XINLINE nrmlp(): lp(2) {} + MSHADOW_XINLINE nrmlp(double l): lp(l) {} + + /* \brief power for Lp norm */ + MSHADOW_XINLINE static double lp_power(volatile double src, volatile double p) { + if (p != 0.0) { + if (src == 0.0) { + return src; + } else { + return power::Map(src, p); + } + } else { // 0-norm, sparsity + return static_cast<double>(src != 0); + } + } + + /*! \brief do reduction into dst */ + template<typename AType, typename DType> + MSHADOW_XINLINE void Reduce(volatile AType& sum_of_powers, volatile DType src) { // NOLINT(*) + if (src != 0) { + sum_of_powers += AType(lp_power(static_cast<double>(src), lp)); + } + } + + /*! \brief do stable reduction into dst */ + template<typename AType, typename DType> + MSHADOW_XINLINE void Reduce(volatile AType& sum_of_powers, volatile DType src, volatile DType& scale) { // NOLINT(*) + if (src != 0) { + DType abs = abs::Map(src); + if (scale < abs) { + sum_of_powers = 1 + sum_of_powers * AType(lp_power(static_cast<double>(scale / abs), lp)); + scale = abs; + } else { + sum_of_powers = sum_of_powers + AType(lp_power(static_cast<double>(abs / scale), lp)); + } + } + } + + /*! \brief combine the results of two reducers */ + template<typename DType> + MSHADOW_XINLINE static void Merge(volatile DType& dst_val, volatile DType& src_val) { // NOLINT(*) + dst_val += src_val; + } + + /*! \brief combine the results of two reducers */ + template<typename DType> + MSHADOW_XINLINE static void Merge(volatile DType& dst_ssq, volatile DType& dst_scale, volatile DType& src_ssq, volatile DType& src_scale) { // NOLINT(*) + if (dst_scale != 0 && dst_scale >= src_scale) { + dst_ssq = dst_ssq + src_ssq * DType(lp_power(static_cast<double>(src_scale / dst_scale), 2)); + } else if (src_scale != 0 && dst_scale < src_scale) { + dst_ssq = src_ssq + dst_ssq * DType(lp_power(static_cast<double>(dst_scale / src_scale), 2)); + dst_scale = src_scale; + } + } + + /*! \brief finalize reduction result */ + template<typename DType> + MSHADOW_XINLINE void Finalize(volatile DType& sum_of_powers) { // NOLINT(*) + if (lp != 0.0) { + sum_of_powers = DType(lp_power(static_cast<double>(sum_of_powers), 1.0 / lp)); + } + } + + /*! \brief finalize reduction result */ + template<typename DType> + MSHADOW_XINLINE void Finalize(volatile DType& sum_of_powers, volatile DType& scale) { // NOLINT(*) + if (lp != 0.0) { + sum_of_powers = scale * DType(lp_power(static_cast<double>(sum_of_powers), 1.0 / lp)); + } + } + + /*! + *\brief set the initial value during reduction + */ + template<typename DType> + MSHADOW_XINLINE static void SetInitValue(DType &sum_of_powers) { // NOLINT(*) + sum_of_powers = 0; + } + + /*! + *\brief set the initial value during reduction + */ + template<typename DType> + MSHADOW_XINLINE static void SetInitValue(DType &sum_of_powers, DType &scale) { // NOLINT(*) + SetInitValue(sum_of_powers); + scale = 0; + } +}; + +/*! \brief Elementwise gradient of Lp-norm, does not handle p = 1 */ +struct nrmlp_grad : public mxnet_op::tunable { + double lp; + MSHADOW_XINLINE nrmlp_grad(): lp(2) {} + MSHADOW_XINLINE nrmlp_grad(double l): lp(l) {} + + /* \brief elementwise gradient of lp norm */ + template<typename DType> + MSHADOW_XINLINE DType Map(DType a, DType b) { + DType ret; + if (lp != 0.0) { // dx_i = (|x_i| / y) ^ (p - 1) * sgn(x_i) + DType abs_a = DType(abs::Map(a)); + DType sgn_a = DType(sign::Map(a)); + ret = power::Map(DType(abs_a / b), DType(lp - 1)) * sgn_a; + } else { // L0 norm is elementwise constant and non-differentiable + ret = 0; + } + return ret; + } +}; + +/*! \brief Gradient for abs-min/max */ +struct abs_grad : public mxnet_op::tunable { + template<typename DType> + MSHADOW_XINLINE static DType Map(DType a, DType b) { + DType sgn = DType(sign::Map(a)); + DType grad = DType(abs::Map(a)) == DType(abs::Map(b)) ? + DType(1.0) : DType(0.0); + return sgn * grad; + } +}; + +/*! \brief Sign */ +struct abs_sign : public mxnet_op::tunable { + template<typename DType> + MSHADOW_XINLINE static DType Map(DType a, DType b) { + return DType(sign::Map(a)); + } +}; + +} // namespace mshadow_op + +inline bool NumpyLpNormShape(const nnvm::NodeAttrs& attrs, + mxnet::ShapeVector *in_attrs, + mxnet::ShapeVector *out_attrs); + +inline bool NumpyMatrixNormShape(const nnvm::NodeAttrs& attrs, + mxnet::ShapeVector *in_attrs, + mxnet::ShapeVector *out_attrs); + +inline void assign_svd_empty(mxnet::ShapeVector *out_attrs); + +bool NumpyNormShape(const nnvm::NodeAttrs& attrs, + mxnet::ShapeVector *in_attrs, + mxnet::ShapeVector *out_attrs); + +TShape swapMatDims(const TShape &shape, const TShape &axis); + +TShape inverseTranspose(const TShape &axes); + +struct NumpyNormParam : public dmlc::Parameter<NumpyNormParam> { + double ord; + dmlc::optional<mxnet::TShape> axis; + bool keepdims; + int flag; + DMLC_DECLARE_PARAMETER(NumpyNormParam) { + DMLC_DECLARE_FIELD(ord).set_default(2) + .describe("Order of the norm. inf means numpy’s inf object."); + DMLC_DECLARE_FIELD(axis).set_default(dmlc::optional<mxnet::TShape>()) + .describe(R"code(If axis is an integer, it specifies the axis of x along + which to compute the vector norms. If axis is a 2-tuple, + it specifies the axes that hold 2-D matrices, and the matrix norms of + these matrices are computed. If axis is None then either a vector norm (when x is 1-D) + or a matrix norm (when x is 2-D) is returned.If axis is an integer, + it specifies the axis of x along which to compute the vector norms. + If axis is a 2-tuple, it specifies the axes that hold 2-D matrices, + and the matrix norms of these matrices are computed. If axis is None then either a + vector norm (when x is 1-D) or a matrix norm (when x is 2-D) is returned.)code"); + DMLC_DECLARE_FIELD(keepdims).set_default(false) + .describe("If this is set to `True`, the reduced axis is left " + "in the result as dimension with size one."); + DMLC_DECLARE_FIELD(flag).set_default(-1) + .describe("Mapping relations between ord and flag." + "ord: None, 'fro', 'nuc', 'inf' '-inf'." + "flag: 0 , 1, 2, 3, 4. "); + } +}; + +TShape swapMatDims(const TShape &shape, const TShape &axis); +TShape inverseTranspose(const TShape &axes); + +bool NumpyNormShape(const nnvm::NodeAttrs& attrs, + mxnet::ShapeVector *in_attrs, + mxnet::ShapeVector *out_attrs); + +template<typename xpu> +void NumpyLpNormCompute(const nnvm::NodeAttrs& attrs, + const OpContext& ctx, + const std::vector<TBlob>& inputs, + const std::vector<OpReqType>& req, + const std::vector<TBlob>& outputs) { + const NumpyNormParam& param = nnvm::get<NumpyNormParam>(attrs.parsed); + double ord = param.ord; + + if (req[0] == kNullOp) return; + + mxnet::TShape small; + mxnet::TShape out_shape = outputs[0].shape_; + if (param.keepdims) { + small = outputs[0].shape_; + } else { + small = ReduceAxesShapeImpl(inputs[0].shape_, param.axis, true, false); + const_cast<std::vector<TBlob>&>(outputs)[0] = outputs[0].reshape(small); + } + bool safe_acc = dmlc::GetEnv("MXNET_SAFE_ACCUMULATION", false); + if (!safe_acc && inputs[0].type_flag_ == mshadow::kFloat16) { + common::LogOnce("MXNET_SAFE_ACCUMULATION=1 is recommended for LpNorm with float16 inputs. " + "See https://mxnet.apache.org/api/faq/env_var " + "for more details."); + } + if (param.axis.value().ndim() != 2) { // elementwise Lp-norm + if (ord == -std::numeric_limits<double>::infinity()) { // -inf norm + LOG(FATAL) << "Under construction. "; Review comment: The case is actually handled in front-end, these are placeholders in case people want me to move the implementation to the back-end, and would not be reached currently. I will change the message to something like "ord = inf handled in front-end" instead. ---------------------------------------------------------------- This is an automated message from the Apache Git Service. To respond to the message, please log on to GitHub and use the URL above to go to the specific comment. For queries about this service, please contact Infrastructure at: us...@infra.apache.org With regards, Apache Git Services
