bfgs.py 8.6 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
# Copyright (c) 2022 PaddlePaddle Authors. All Rights Reserved.
#
# Licensed 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.

import numpy as np

from .line_search import strong_wolfe
18 19 20 21 22
from .utils import (
    _value_and_gradient,
    check_input_type,
    check_initial_inverse_hessian_estimate,
)
23 24 25 26

import paddle


27 28 29 30 31 32 33 34 35 36 37 38 39
def minimize_bfgs(
    objective_func,
    initial_position,
    max_iters=50,
    tolerance_grad=1e-7,
    tolerance_change=1e-9,
    initial_inverse_hessian_estimate=None,
    line_search_fn='strong_wolfe',
    max_line_search_iters=50,
    initial_step_length=1.0,
    dtype='float32',
    name=None,
):
40 41
    r"""
    Minimizes a differentiable function `func` using the BFGS method.
S
Sing_chan 已提交
42 43 44
    The BFGS is a quasi-Newton method for solving an unconstrained optimization problem over a differentiable function.
    Closely related is the Newton method for minimization. Consider the iterate update formula:

45
    .. math::
S
Sing_chan 已提交
46 47 48
        x_{k+1} = x_{k} + H_k \nabla{f_k}

    If :math:`H_k` is the inverse Hessian of :math:`f` at :math:`x_k`, then it's the Newton method.
49
    If :math:`H_k` is symmetric and positive definite, used as an approximation of the inverse Hessian, then
50
    it's a quasi-Newton. In practice, the approximated Hessians are obtained
51
    by only using the gradients, over either whole or part of the search
S
Sing_chan 已提交
52
    history, the former is BFGS, the latter is L-BFGS.
53

54
    Reference:
S
Sing_chan 已提交
55
        Jorge Nocedal, Stephen J. Wright, Numerical Optimization, Second Edition, 2006. pp140: Algorithm 6.1 (BFGS Method).
56 57

    Args:
58
        objective_func: the objective function to minimize. ``objective_func`` accepts a 1D Tensor and returns a scalar.
59
        initial_position (Tensor): the starting point of the iterates, has the same shape with the input of ``objective_func`` .
S
Sing_chan 已提交
60 61 62
        max_iters (int, optional): the maximum number of minimization iterations. Default value: 50.
        tolerance_grad (float, optional): terminates if the gradient norm is smaller than this. Currently gradient norm uses inf norm. Default value: 1e-7.
        tolerance_change (float, optional): terminates if the change of function value/position/parameter between two iterations is smaller than this value. Default value: 1e-9.
63
        initial_inverse_hessian_estimate (Tensor, optional): the initial inverse hessian approximation at initial_position. It must be symmetric and positive definite. If not given, will use an identity matrix of order N, which is size of ``initial_position`` . Default value: None.
S
Sing_chan 已提交
64 65 66
        line_search_fn (str, optional): indicate which line search method to use, only support 'strong wolfe' right now. May support 'Hager Zhang' in the futrue. Default value: 'strong wolfe'.
        max_line_search_iters (int, optional): the maximum number of line search iterations. Default value: 50.
        initial_step_length (float, optional): step length used in first iteration of line search. different initial_step_length may cause different optimal result. For methods like Newton and quasi-Newton the initial trial step length should always be 1.0. Default value: 1.0.
67
        dtype ('float32' | 'float64', optional): data type used in the algorithm, the data type of the input parameter must be consistent with the dtype. Default value: 'float32'.
S
Sing_chan 已提交
68 69
        name (str, optional): Name for the operation. For more information, please refer to :ref:`api_guide_Name`. Default value: None.

70
    Returns:
S
Sing_chan 已提交
71 72 73 74 75 76 77 78
        output(tuple):

            - is_converge (bool): Indicates whether found the minimum within tolerance.
            - num_func_calls (int): number of objective function called.
            - position (Tensor): the position of the last iteration. If the search converged, this value is the argmin of the objective function regrading to the initial position.
            - objective_value (Tensor): objective function value at the `position`.
            - objective_gradient (Tensor): objective function gradient at the `position`.
            - inverse_hessian_estimate (Tensor): the estimate of inverse hessian at the `position`.
79 80 81 82 83

    Examples:
        .. code-block:: python

            import paddle
84

85 86 87 88
            def func(x):
                return paddle.dot(x, x)

            x0 = paddle.to_tensor([1.3, 2.7])
89
            results = paddle.incubate.optimizer.functional.minimize_bfgs(func, x0)
90 91 92 93 94 95 96 97 98 99
            print("is_converge: ", results[0])
            print("the minimum of func is: ", results[2])
            # is_converge:  is_converge:  Tensor(shape=[1], dtype=bool, place=Place(gpu:0), stop_gradient=True,
            #        [True])
            # the minimum of func is:  Tensor(shape=[2], dtype=float32, place=Place(gpu:0), stop_gradient=True,
            #        [0., 0.])
    """

    if dtype not in ['float32', 'float64']:
        raise ValueError(
100 101 102 103
            "The dtype must be 'float32' or 'float64', but the specified is {}.".format(
                dtype
            )
        )
104 105 106 107 108 109 110 111

    op_name = 'minimize_bfgs'
    check_input_type(initial_position, 'initial_position', op_name)

    I = paddle.eye(initial_position.shape[0], dtype=dtype)
    if initial_inverse_hessian_estimate is None:
        initial_inverse_hessian_estimate = I
    else:
112 113 114 115 116
        check_input_type(
            initial_inverse_hessian_estimate,
            'initial_inverse_hessian_estimate',
            op_name,
        )
117 118 119
        check_initial_inverse_hessian_estimate(initial_inverse_hessian_estimate)

    Hk = paddle.assign(initial_inverse_hessian_estimate)
120 121
    # use detach and assign to create new tensor rather than =, or xk will share memory and grad with initial_position
    xk = paddle.assign(initial_position.detach())
122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144

    value, g1 = _value_and_gradient(objective_func, xk)
    num_func_calls = paddle.full(shape=[1], fill_value=1, dtype='int64')

    # when the dim of x is 1000, it needs more than 30 iters to get all element converge to minimum.
    k = paddle.full(shape=[1], fill_value=0, dtype='int64')
    done = paddle.full(shape=[1], fill_value=False, dtype='bool')
    is_converge = paddle.full(shape=[1], fill_value=False, dtype='bool')

    def cond(k, done, is_converge, num_func_calls, xk, value, g1, Hk):
        return (k < max_iters) & ~done

    def body(k, done, is_converge, num_func_calls, xk, value, g1, Hk):
        #############    compute pk    #############
        pk = -paddle.matmul(Hk, g1)

        #############    compute alpha by line serach    #############
        if line_search_fn == 'strong_wolfe':
            alpha, value, g2, ls_func_calls = strong_wolfe(
                f=objective_func,
                xk=xk,
                pk=pk,
                initial_step_length=initial_step_length,
145 146
                dtype=dtype,
            )
147 148
        else:
            raise NotImplementedError(
149 150 151 152
                "Currently only support line_search_fn = 'strong_wolfe', but the specified is '{}'".format(
                    line_search_fn
                )
            )
153 154 155 156 157 158 159 160 161 162 163 164 165 166
        num_func_calls += ls_func_calls

        #############    update Hk    #############
        sk = alpha * pk
        yk = g2 - g1

        xk = xk + sk
        g1 = g2

        sk = paddle.unsqueeze(sk, 0)
        yk = paddle.unsqueeze(yk, 0)

        rhok_inv = paddle.dot(yk, sk)
        rhok = paddle.static.nn.cond(
167
            rhok_inv == 0.0,
168
            lambda: paddle.full(shape=[1], fill_value=1000.0, dtype=dtype),
169 170
            lambda: 1.0 / rhok_inv,
        )
171 172 173

        Vk_transpose = I - rhok * sk * yk.t()
        Vk = I - rhok * yk * sk.t()
174 175 176 177
        Hk = (
            paddle.matmul(paddle.matmul(Vk_transpose, Hk), Vk)
            + rhok * sk * sk.t()
        )
178 179 180 181 182 183

        k += 1

        #############    check convergence    #############
        gnorm = paddle.linalg.norm(g1, p=np.inf)
        pk_norm = paddle.linalg.norm(pk, p=np.inf)
184
        paddle.assign(
185 186
            done | (gnorm < tolerance_grad) | (pk_norm < tolerance_change), done
        )
187 188
        paddle.assign(done, is_converge)
        # when alpha=0, there is no chance to get xk change.
189
        paddle.assign(done | (alpha == 0.0), done)
190 191 192 193 194
        return [k, done, is_converge, num_func_calls, xk, value, g1, Hk]

    paddle.static.nn.while_loop(
        cond=cond,
        body=body,
195 196
        loop_vars=[k, done, is_converge, num_func_calls, xk, value, g1, Hk],
    )
197
    return is_converge, num_func_calls, xk, value, g1, Hk