最速下降优化算法的C++实现

最新推荐文章于 2024-08-04 21:27:58 发布
最新推荐文章于 2024-08-04 21:27:58 发布
阅读量337 收藏
点赞数
文章标签： c/c++ 操作系统 python
原文链接：https://my.oschina.net/zmjerry/blog/11025
版权

2019独角兽企业重金招聘Python工程师标准>>> hot3.png

头文件：

/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                                objfunc.h
 *
 * Function Object.
 *
 * We can use this functor computing the value of objective function and it's
 * gradient vector. The objective function is supposed to be multidimensional,
 * one dimention is the special case of "vector.dim()=1".
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


#ifndef OBJFUNC_H
#define OBJFUNC_H


#include <vector.h>


namespace splab
{

    template <typename Type>
    class ObjFunc
    {

    public:

        /**
         * Initialize the parameters
         */
        ObjFunc( Type aa, Type bb, Type cc ) : a(aa), b(bb), c(cc)
        { }

        /**
         * Compute the value of objective function at point x.
         */
        Type operator()( Vector<Type> &x )
        {
            return ( a*x(1) * exp(b*x(1)*x(1)+c*x(2)*x(2)) );
        }

        /**
         * Compute the gradient of objective function at point x.
         */
        Vector<Type> grad( Vector<Type> &x )
        {
            Type expTerm = exp(a*x(1)*x(1)+b*x(2)*x(2));
            Vector<Type> df(x.dim());
            df(1) = (a+2*a*b*x(1)*x(1)) * expTerm;
            df(2) = 2*a*c*x(1)*x(2) * expTerm;

            return df;
        }

    private:

        // parameters
        Type a,
             b,
             c;

    };
    // class ObjFunc

}
// namespace splab


#endif
// OBJFUNC_H
/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                               linesearch.h
 *
 * Line search algorithm based on quadratic polynomial interpolation.
 *
 * This routine is often used in orther optimation method, so it is designed
 * as a templated base class.
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


#ifndef LINESEARCH_H
#define LINESEARCH_H


#include <vector.h>


namespace splab
{

    template <typename Dtype, typename Ftype>
    class LineSearch
    {

     public:

        LineSearch();
        ~LineSearch();

        Dtype getStep( Ftype &func, Vector<Dtype> &xk, Vector<Dtype> &dk,
                       int maxItr=10 );

        int getFuncNum() const;
        bool isSuccess() const;

    protected:

        int     funcNum;
        bool    success;

    };
    // class LineSearch


    #include <linesearch-impl.h>

}
// namespace splab


#endif
// LINESEARCH_H
/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                               steepdesc.h
 *
 * Steepest descent optimal method( gradient algorithm ).
 *
 * This class is designed for finding the minimum value of objective function
 * in one dimension or multidimension. Inexact line search algorithm is used
 * to get the step size in each iteration.
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


#ifndef STEEPDESC_H
#define STEEPDESC_H


#include <linesearch.h>


namespace splab
{

    template <typename Dtype, typename Ftype>
    class SteepDesc : public LineSearch<Dtype, Ftype>
    {

    public:

        SteepDesc();
        ~SteepDesc();

        void optimize( Ftype &func, Vector<Dtype> &x0,
                       Dtype tol=Dtype(1.0e-6), int maxItr=100 );

        Vector<Dtype> getOptValue() const;
        Vector<Dtype> getGradNorm() const;
        Dtype getFuncMin() const;
        int getItrNum() const;

    private:

        // iteration number
        int itrNum;

        // minimum value of objective function
        Dtype fMin;

        // optimal solution
        Vector<Dtype> xOpt;

        // gradient norm for each iteration
        Vector<Dtype> gradNorm;

    };
    // class SteepDesc


    #include <steepdesc-impl.h>

}
// namespace splab


#endif
// STEEPDESC_H

实现文件：

/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                               linesearch-impl.h
 *
 * Implementation for LineSearch class.
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


/**
 * constructors and destructor
 */
template <typename Dtype, typename Ftype>
LineSearch<Dtype, Ftype>::LineSearch()
{
    funcNum = 0;
    success = true;
}

template <typename Dtype, typename Ftype>
LineSearch<Dtype, Ftype>::~LineSearch()
{
}


/**
 * Finding the step size at point "xk" in direction of "dk" of function
 * "func". The default heuristics number is "maxItr=10".
 */
template <typename Dtype, typename Ftype>
Dtype LineSearch<Dtype, Ftype>::getStep( Ftype &func, Vector<Dtype> &xk,
                                         Vector<Dtype> &dk, int maxItr )
{
    // Set line search parameters that everyone uses.
    Dtype mu = Dtype(0.001),
          kUp = Dtype(0.5),
          kLow = Dtype(0.1),
          alpha = Dtype(1.0),
          alphaMin,
          alphaMax;

    Dtype fNew,
          fk = func(xk);

    Vector<Dtype> xNew,
                  gk = func.grad(xk);

    Dtype gd = dotProd( gk, dk );

    for( int i=0; i<maxItr; ++i )
    {
        xNew = xk + alpha*dk;
        fNew = func(xNew);
        funcNum++;

        if( fNew < fk+mu*alpha*gd )
        {
            success = true;
            return alpha;
        }
        else
        {
            alphaMin = kLow*alpha;
            alphaMax = kUp*alpha;

            // Compute the step by using quadratic polynomial interpolation.
            alpha = Dtype(-0.5)*alpha*alpha*gd / ( fNew-fk-alpha*gd );

            // bound checking
            if( alpha < alphaMin )
                alpha = alphaMin;
            else if( alpha > alphaMax )
                alpha = alphaMax;
        }
    }

    if( fNew>=fk )
    {
        success = false;
        return Dtype(0.0);
    }
    else
    {
        success = true;
        return alpha;
    }
}


/**
 * Get the number of objective function's calculation.
 */
template <typename Dtype, typename Ftype>
inline int LineSearch<Dtype, Ftype>::getFuncNum() const
{
    return funcNum;
}


/**
 * Judgement whether the optimal solution is found or not.
 */
template <typename Dtype, typename Ftype>
inline bool LineSearch<Dtype, Ftype>::isSuccess() const
{
    return success;
}
/*
 * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 2 or any later version.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 *
 * 1. Redistributions of source code must retain the above copyright notice,
 *    this list of conditions and the following disclaimer.
 *
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
 * more details. A copy of the GNU General Public License is available at:
 * http://www.fsf.org/licensing/licenses
 */


/*****************************************************************************
 *                               steepdesc-impl.h
 *
 * Implementation for SteepDesc class.
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


/**
 * constructors and destructor
 */
template <typename Dtype, typename Ftype>
SteepDesc<Dtype, Ftype>::SteepDesc() : LineSearch<Dtype, Ftype>()
{
}

template <typename Dtype, typename Ftype>
SteepDesc<Dtype, Ftype>::~SteepDesc()
{
}


/**
 * Finding the optimal solution. The default tolerance error and maximum
 * iteratin number are "tol=1.0e-6" and "maxItr=100", respectively.
 */
template <typename Dtype, typename Ftype>
void SteepDesc<Dtype, Ftype>::optimize( Ftype &func, Vector<Dtype> &x0,
                                        Dtype tol, int maxItr )
{
    // initialize parameters.
    int k = 0;
    Vector<Dtype> x(x0);
    Dtype fx = func(x);
    this->funcNum++;
    Vector<Dtype> gnorm(maxItr);
    Vector<Dtype> g = func.grad(x);
    gnorm[k++]= norm(g);

    Dtype alpha;
    Vector<Dtype> d(g.dim());

    while( ( gnorm(k) > tol ) && ( k < maxItr ) )
    {
        // descent direction
        d =  - g;

        // one dimension searching
        alpha = this->getStep( func, x, d );
        if( !this->success )
            break;

		// update
        x += alpha*d;
        fx = func(x);
        this->funcNum++;
        g = func.grad(x);
        gnorm[k++] = norm(g);
    }

    xOpt = x;
    fMin = fx;
    gradNorm.resize(k);
    for( int i=0; i<k; ++i )
        gradNorm[i] = gnorm[i];

    if( gradNorm[k-1] > tol )
        this->success = false;
}


/**
 * Get the optimum point.
 */
template <typename Dtype, typename Ftype>
inline Vector<Dtype> SteepDesc<Dtype, Ftype>::getOptValue() const
{
    return xOpt;
}


/**
 * Get the norm of gradient in each iteration.
 */
template <typename Dtype, typename Ftype>
inline Vector<Dtype> SteepDesc<Dtype, Ftype>::getGradNorm() const
{
    return gradNorm;
}


/**
 * Get the minimum value of objective function.
 */
template <typename Dtype, typename Ftype>
inline Dtype SteepDesc<Dtype, Ftype>::getFuncMin() const
{
    return fMin;
}


/**
 * Get the iteration number.
 */
template <typename Dtype, typename Ftype>
inline int SteepDesc<Dtype, Ftype>::getItrNum() const
{
    return gradNorm.dim()-1;
}

测试代码：

/*****************************************************************************
 *                               steepdesc_test.cpp
 *
 * Steepest descent method testing.
 *
 * Zhang Ming, 2010-03, Xi'an Jiaotong University.
 *****************************************************************************/


#define BOUNDS_CHECK

#include <iostream>
#include <iomanip>
#include <objfunc.h>
#include <steepdesc.h>


using namespace std;
using namespace splab;


typedef float  Type;


int main()
{
    Type a = 1.0,
         b = -1.0,
         c = -1.0;
    ObjFunc<Type> f( a, b, c );
    Vector<Type> x0(2);
    x0(1) = Type(0.0);
    x0(2) = Type(0.0);

    Type tolErr = 1.0e-3;
    SteepDesc< Type, ObjFunc<Type> > steep;
    steep.optimize( f, x0, tolErr );
    if( steep.isSuccess() )
    {
        Vector<Type> xmin = steep.getOptValue();
        int N = steep.getItrNum();
        cout << "The iterative number is:   " << N << endl << endl;
        cout << "The number of function calculation is:   "
             << steep.getFuncNum() << endl << endl;
        cout << setiosflags(ios::fixed) << setprecision(4);
        cout << "The optimal value of x is:   " << xmin << endl;
        cout << "The minimum value of f(x) is:   " << f(xmin) << endl << endl;
        cout << "The gradient's norm at x is:   "
             << steep.getGradNorm()[N] << endl << endl;
    }
    else
        cout << "The optimal solution  can't be found!" << endl;

    return 0;
}

运行结果：

The iterative number is:   14

The number of function calculation is:   43

The optimal value of x is:   size: 2 by 1
-0.7070
0.0000

The minimum value of f(x) is:   -0.4289

The gradient's norm at x is:   0.0006


Process returned 0 (0x0)   execution time : 0.062 s
Press any key to continue.

转载于:https://my.oschina.net/zmjerry/blog/11025

weixin_33953384
关注
  • 0
    点赞
  • 踩
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
C语言——自定义函数向右移动数组中的元素
qq_59990833的博客
08-01 343
C语言——自定义函数向右移动数组中的元素
【C语言】Top K问题【建小堆】
2301_80840905的博客
08-04 300
在我们生活中，经常会遇到TopK问题。建小堆解决
C++ STL copy, move 用法
08-04 132
前向（从前向后的顺序）拷贝/移动操作，将一个容器元素拷贝/移动到另一容器中。
C++ | string
2302_78381559的博客
07-31 625
在C++中，类是用于处理文本字符串的一个非常强大的工具。它定义在头文件中，并且是标准模板库 (STL) 的一部分。提供了一种方便的方式来存储和操作字符序列，使得开发人员能够轻松地执行常见的字符串操作，而不需要直接管理内存。auto在 C++11 中用于简化类型声明，使代码更简洁。使用auto时要注意变量类型的一致性和正确性。auto不能用于声明函数参数或直接声明数组。
C语言指针详解（四）
Duan1044592637的博客
08-01 768
1、冒泡排序2、回调函数（举例通过转移表实现简易计算器）3、qsort函数使用举例4、采用冒泡排序方式进行qsort函数模拟实现
【doghead】h264测试文件、读取、模拟采集时间戳及packtizer
最新发布
突围
08-04 26
dog head 模拟h264的采集
C++ 继承 派生类的拷贝构造
gopher9511的博客
08-04 448
拷贝构造，也是一种构造函数，也没有被继承下来。
C语言程序设计之链表篇2
繁华落尽心殇
08-02 905
建立单向链表的一般步骤是：建立头指针 -> 建立第一个结点 -> 头指针指向第一个结点 -> 建立第二个结点 -> 头指针指向第二个结点。考查链表的数据结构，需利用指针变量才能实现，一个结点中应包含一个指针变量，用来存放下一个结点的地址。的功能是：将不带头结点的单向链表结点数据中的数据从小到大排序。的功能是：将不带头结点的单向链表逆置，即若原链表中从头自尾结点数据域依次为。的功能是 ：统计带头结点的单向链表中结点的个数，并存放在形参。在带头结点的单向链表中，查找数据域中值为。的结点，函数返回 0 值。
「C++系列」字符串
日常笔记/总结/系列知识梳理
08-01 963
在C++中，字符串是一种用于存储和操作文本数据的类型。C++标准库提供了几种处理字符串的方式，但最常用的是类，它定义在头文件中。提供了丰富的成员函数来执行各种字符串操作，比如插入、删除、查找、替换、连接字符串等。
[C++] 模板编程-10 模板形参----＞模板参数
qq_45280097的博客
07-30 214
首先我们里面的template相当于一个模板, 后面的class与typename功能一样,表示指定类型名为T，而这里的类型是前面的template
C++理解“抛出一个异常”与“传递一个参数”或“调用一个虚函数”间的差异
changjian198384的博客
08-02 544
不论通过传值捕获异常还是通过引用捕获（不能通过指针捕获这个异常，因为类型不匹配）都将进行lcalWidget的拷贝操作，也就说传递到catch子句中的是localWidget的拷贝。如果把localWidget本身（而不是它的拷贝）传递给catch子句，这个子句接收到的只是一个被析构了的Widget，一个Widget的“尸体”。（因为这样的代码在C++里通常是不合法的。这两个catch块的差别在于第一个catch块中重新抛出的是当前捕获的异常，而第二个catch块中重新抛出的是当前捕获异常的一个新的拷贝。
深入理解 C++ 中的 const 成员函数
威桑的博客
08-03 313
const 成员函数
c语言-结构体
m0_71813740的博客
07-31 674
a.普通的定义b.宏定义-- #define-- 它是预处理中的宏定义命令，在预处理时进行简单的字符串替换，不作正确性检测。-- #define的弊端：-- 只有a被定义为指针类型，b是int类型。这是因为宏定义只是简单的替换，很容易再使用的时候出现错误，如放在式子中优先级的问题之类。c.类型重定义-- typedef-- typedef是关键字，在编译时处理，有类型检测的功能，它在自己的作用域内给一个数据类型取一个新的名字。typedef主要用来定义数组、指针、结构体等类型。2、初始化。
对于n进制转十进制的解法及代码（干货！）
yykYYK2010的博客
08-01 157
对于p进制转十进制，我们有：(x)p=a[0]*p^0+a[1]*p^1+a[2]*p^2+...+a[n]*p^n。据此，我们可以编出c++代码来解决问题。举个例子：（11001）
【LeetCode Cookbook（C++ 描述）】一刷哈希表（Hash Table）（下）
David Max 的博客
08-02 882
本系列文章仅是 GitHub 大神 @halfrost 的刷题笔记《LeetCode Cookbook》的提纲以及示例、题集的 C++转化。原书请自行下载学习。本篇文章涉及新手应该优先刷的几道经典哈希表算法题，分两篇文章，以后会更新“二刷”“三刷”等等。
C++ 基础练习 - Chapter 8 （英文版）
weixin_35193147的博客
07-30 702
多多练习，做多思考。加油！
八：C语言-函数的基本概念
qiao_yue的博客
08-01 1056
本篇文章详细讲解了库函数与自定义函数的区别、函数的嵌套调用及链式访问、函数的声明和定义、static和extern等基础知识
C++为何难以割舍历史包袱？
T1246632660的博客
08-01 413
你以为的历史包袱实际上是C++几十年以来的丰厚遗产，而且其中有的遗产非常有价值，还有有相当多的代码非常稳定，如果使用这些代码的人不是追新强迫真，它们永远都不需要重写。C++已经占领了各种工业软件的几乎每个角落，这其中牵涉到每个行业里面大量的数据结构、算法、以及为海量的论文所实现的代码，这些代码完全重写所需要大量的时间和金钱，而且现在的人力成本远超当初，所以完全重写几乎已经不可能。C++的代码到底有多少行，这个没人清楚，因为很多不开源。BS在演讲的时候说过是bilions。兼容对普通人有啥意义呢?
C++ 学习(2) ---- std::cout 格式化输出
小猪佩奇的博客
08-04 483
C++ 通常使用cout输出数据，和printf()函数相比，cout实现格式化输出数据的方式更加多样化；一方面，cout 作为ostream类的对象，该类中提供有一些成员方法，可实现对输出数据的格式化；另一方面，为了方面用户格式化输出数据，C++ 标准库专门提供了一个 <iomanip> 头文件，该头文件中包含有大量的格式控制符，也就是流操作算子，使用上更加方便；
KMP算法详解：C++实现与性能优化
C++实现KMP算法时，首先需要初始化next数组。这通常通过遍历模式串，逐个比较字符并更新next数组完成。接着，在匹配过程中，通过比较文本串和模式串的字符，如果匹配则继续，如果不匹配，则根据next数组的值调整模式...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值