查看: 5151|回复: 3|关注: 0
应用sym函数求方程解析解的问题(新的问题)
本帖最后由 stilljourney 于
05:58 编辑
问题描述：u是[0, 1.2]之间的自变量，p是u的函数，其关系可以通过设一个矩阵的行列式为零来表示，即det[D]=0，这个式子具有这样的a*p^4+b*p^3+b*p^2+d*p=0的多项式形式，p有四个解，分别可以用P=Pr+i*Pi表示。
目的：得到四个p值关于u的表达式，在u=[0, 1.2]的范围内分别画出其实部、虚部关于u的曲线
我的思路：用sym函数写出来det[D]=0的式子，用[p]=solve(equation, u)的命令求解p关于u的解析解表达式
我目前写出的程序：
function FlutterAnalysis(dampingOrNot)
%Set numerical values to parameters
c=2; %chord
b=c/2; %length scale (half-chord)
M=1; %mass
S_alpha=0.05; %static unbalance
I_alpha=0.25; %moment of inertia
K_alpha=0.25; %torsional stiffness
omega_alpha=K_alpha/I_
K_h=0.25; %translational stiffness
omega_h=K_h/M;
rho=1/(10*pi); %density of the incoming flow
aeroCoef=2* %aerodynamic coefficient
e=b*0.4; %the distance between the aerodynamic center and the elastic axis
switch dampingOrNot
& & case 1
& && &&&%Create det[D]=0 function to calculate p
& && &&&equation=(M*(p^2)+K_h)*(I_alpha*(p^2)+K_alpha-e*0.5*rho*(u^2)* ...
& && && && &c*aeroCoef) - (S_alpha*(p^2)+0.5*rho*(u^2)*c*aeroCoef)* ...
& && && && &S_alpha*(p^2);
& && &&&[p]=solve(equation, p);
& & case 2
& && &&&%Create det[D]=0 function to calculate p
& && &&&equation=((M*(p^2)+K_h+0.5*rho*(u^2)*c*aeroCoef*(p/u))* ...
& && && && &(I_alpha*(p^2)+K_alpha-e*0.5*rho*(u^2)*c*aeroCoef))- ...
& && && && &((S_alpha*(p^2)+0.5*rho*(u^2)*c*aeroCoef)* ...
& && && && &(S_alpha*(p^2)-e*0.5*rho*(u^2)*c*aeroCoef*(p/u)));
& && &&&[p]=solve(equation, p)
& & xaxis=0:0.1:1.2;
& & u=xaxis*omega_
& & yaxis_real=real(subs(p(i, :)));
& & yaxis_imag=imag(subs(p(i, :)));
& & axis([0, 1.5, -1.2, 1.2]);
& & plot(xaxis,yaxis_real, '-', 'LineWidth', 2);
& & plot(xaxis,yaxis_imag, '-', 'LineWidth', 2);
& & hold all
目前的问题：case 1可以计算成功，但是case 2就会出现问题，可能是因为方程变复杂了哪里出错了，运行case 2的结果是：
&& FlutterAnalysis(2)
RootOf(z^4 + (4065*pi*u*z^3)/68928 - (400*z^2*((03147*pi*u^2)/8558720 - 5/16))/99 + (453175*pi*u*z)/8716416 - (3475*pi*u^2)/301824 + 25/99, z)
??? Error using ==& plot
Conversion to double from sym is not possible.
Error in ==& FlutterAnalysis at 45
& & plot(xaxis,yaxis_real, '-', 'LineWidth', 2);
不知道z是哪里跑出来的，图像也是一片空白，并且还有一行报错，在运行case 1之后就不会有报错。
继续求助，谢谢！！
关注者： 82
RootOf (f(z), z)的意思是说，“p的值恰好是这个关于z的多项式方程的根”。你可以理解为p就是z。
MATLAB默认只解三次或以下的方程，四次以上的方程根都是用RootOf的形式给出的。你可以在solve里设置MaxDegree为4，这样就可以解出这个四次方程的解析解。关于如何设置，请查阅官方文档。我设置MaxDegree没有成功过，不知道为什么。
stellari 发表于
RootOf (f(z), z)的意思是说，“p的值恰好是这个关于z的多项式方程的根”。你可以理解为p就是z。
MATLAB默 ...
谢谢stellari！我试了，也查了mathwork文档，确实也没有解成功。似乎这条路就有点走不通了
对这样一个多项式求值问题，你有没有什么建议能尝试些别的什么方法呢？
关注者： 1
你的程序让我觉得很复杂。我的matlab怎么学不好啊
站长推荐 /1
Powered bysyms xsolve(sqrt(1/x+1)+sqrt(1/(x-1))-sqrt(1/(x^2-2))-sqrt(x^2-1))这个方程是存在根式解的(应该是5个解)，但是以上代码却只返回了一个数值解。请问，怎么得到解析解呢？如果确实得不到解析解，怎么全局求解所有的数值解(fzero只能局部求解)？
用Mathematica用matlab求解的话，可以定义不同的初始值进行尝试，比如-1，1，i，-i等
数值求解就这样
用Mathematica做的，偷懒直接一图流了，不喜勿怪&br&&img src=&/d1afbfb47c_b.png& data-rawwidth=&662& data-rawheight=&1080& class=&origin_image zh-lightbox-thumb& width=&662& data-original=&/d1afbfb47c_r.png&&&br&===================更新======================&br&@&a href=&/people/1235813& class=&internal&&曹洪洋&/a&指出用Reduce或者Solve加Method-&Reduce能很快算出所有的解析解，效果如下&br&&img src=&/71a2f04eea05d_b.png& data-rawwidth=&1628& data-rawheight=&208& class=&origin_image zh-lightbox-thumb& width=&1628& data-original=&/71a2f04eea05d_r.png&&用时不足一秒，果然不知高到哪里去了
用Mathematica做的，偷懒直接一图流了，不喜勿怪===================更新======================@指出用Reduce或者Solve加Method-&Reduce能很快算出所有的解析解，效果如下用时不足一秒，果然不知高到哪里去了
这个方程初步可以写成&img src=&///equation?tex=x%5E2-2-%5Cfrac%7B1%7D%7Bx%7D%2B%5Cfrac%7B1%7D%7Bx-1%7D-%5Cfrac%7B1%7D%7Bx%5E2-2%7D%3D2%5Csqrt%7Bx%2B1%7D-2%5Csqrt%7B%5Cfrac%7Bx%2B1%7D%7Bx%7D%5Cfrac%7B1%7D%7Bx%5E2-2%7D%7D& alt=&x^2-2-\frac{1}{x}+\frac{1}{x-1}-\frac{1}{x^2-2}=2\sqrt{x+1}-2\sqrt{\frac{x+1}{x}\frac{1}{x^2-2}}& eeimg=&1&&，目测至少超过了5次，理论上是不能一般地用根式表达出来的，机器更做不到这一点了。solve函数只在3次以下可以有解析解。不过要求数值解的话可以用fsolve函数，这个函数可以先赋初值来得到最优解，比如初值赋为-2的解是这样的&br&&div class=&highlight&&&pre&&code class=&language-text&&x = fsolve(@(x) sqrt(1/x+1)+sqrt(1/(x-1))-sqrt(1/(x^2-2))-sqrt(x^2-1),[-2],optimoptions('fsolve'))
&/code&&/pre&&/div&一般情况下不论赋多大的初值，算出始终是同样的解。不过这个方程有若干个间断点，赋某个初值可能只能够求到一段上的解，因为当计算机运行到断点上时会自动停止运算。所以要找出其间断点，然后分开来求。&br&在这个式子上可以看到其断点分别是在&img src=&///equation?tex=x%3D0%2Cx%3D1%2Cx%3D%5Csqrt%7B2%7D& alt=&x=0,x=1,x=\sqrt{2}& eeimg=&1&&时。所以只需在小于&img src=&///equation?tex=0& alt=&0& eeimg=&1&&,&img src=&///equation?tex=0& alt=&0& eeimg=&1&&到&img src=&///equation?tex=1& alt=&1& eeimg=&1&&，大于&img src=&///equation?tex=1& alt=&1& eeimg=&1&&，&img src=&///equation?tex=1& alt=&1& eeimg=&1&&到&img src=&///equation?tex=%5Csqrt%7B2%7D& alt=&\sqrt{2}& eeimg=&1&&这几个范围内各取几个有理数作为初值来算就行了。另外&img src=&///equation?tex=-2& alt=&-2& eeimg=&1&&到&img src=&///equation?tex=0& alt=&0& eeimg=&1&&之间解可能比较密集，函数对初值比较敏感，所以要多取几次，最好写一个循环。
这个方程初步可以写成x^2-2-\frac{1}{x}+\frac{1}{x-1}-\frac{1}{x^2-2}=2\sqrt{x+1}-2\sqrt{\frac{x+1}{x}\frac{1}{x^2-2}}，目测至少超过了5次，理论上是不能一般地用根式表达出来的，机器更做不到这一点了。solve函数只在3次以下可以有解析解。不过要求…
已有帐号？
无法登录？
社交帐号登录查看: 2252|回复: 1
MATLAB软件解密——绝好的数据分析软件使用
在线时间912 小时
阅读权限60
该用户从未签到
主题帖子积分
银牌会员, 积分 2242, 距离下一级还需 158 积分
本帖最后由 xlong87 于
15:58 编辑
MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
基本功能　　MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、、科学数据可视化以及非动态系统的和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。　
　　MATLAB和、并称为三大软件。它在数学类科技应用软件中在方面首屈一指。MATLAB可以进行运算、绘制和数据、实现、创建用户界面、连 、、金融建模设计与分析等领域。
　　MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的。在新的版本中也加入了对，，，的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。
基本应用　　MATLAB 产品族可以用来进行以下各种工作：
　　●数值分析
　　●数值和符号计算
　　●工程与科学绘图
　　●控制系统的设计与仿真
　　●技术
　　●技术
　　●通讯系统设计与仿真
MATLAB在通讯系统设计与仿真的应用
　　●财务与金融工程
　　●管理与调度优化计算（运筹学）
　　MATLAB 的应用范围非常广，包括信号和图像处理、通讯、控制系统设计、测试和测量、财务建模和分析以及计算学等众多应用领域。附加的工具箱（单独提供的专用MATLAB 函数集）扩展了MATLAB 环境，以解决这些应用领域内特定类型的问题。
13:40 上传
点击文件名下载附件
下载积分: 粮票 -1
7.23 MB, 下载次数: 9, 下载积分: 粮票 -1
途之人可以为之禹
在线时间17 小时
阅读权限20
该用户从未签到
主题帖子积分
谢啦。。。
杰出参与勋章
杰出参与勋章
Powered by评论-2096&
trackbacks-0
　　本节主要是练习下斯坦福DL网络教程UFLDL关于Sparse coding那一部分，具体的网页教程参考：中分别采用sparse coding和拓扑的sparse coding方法进行学习，并观察学习到的这些图像基图像的特征。训练数据时自然图片IMAGE，在给出的教程网站上有。
　　实验基础
　　Sparse coding的主要是思想学习输入数据集&基数据&，一旦获得这些&基数据&，输入数据集中的每个数据都可以用这些&基数据&的线性组合表示，而稀疏性则体现在这些线性组合系数是系数的，即大部分的值都为0。很显然，这些&基数据&的尺寸和原始输入数据的尺寸是相同的，另外&基数据&的个数通常要比每个样本的维数大。最简单的理解可以看前面博文提到过的公式：
　　其中的输入数据x可以分解成基Ф的线性组合，ai为组合系数。不过那只是针对一个数据而已，而在ML领域中都是大数据，因此下面来考虑样本是矩阵的形式。在前面博文中我们已经介绍过sparse coding系统非拓扑时的代价函数为：
　　拓扑结构时的代价函数如下：
　　在训练阶段我们的目的是要通过优化算法求出最佳的参数A，因为A就是我们的&基数据&集。但是以上2个代价函数表达式中都有两个未知的参数矩阵，即A和s，所以不能采用简单的优化方法。此时一般的优化思想为交叉优化，即先固定一个A来优化s，然后固定该s来优化A，以此类推，等迭代步骤到达预设值时就停止。而在优化过程中首先要解决的就是代价函数对参数矩阵A和s的求导问题。
　　此时的求导涉及到了矩阵范数的求导，一般有2种方法，第一种是将求导问题转换到矩阵的迹的求导，可以参考前面博文。第二种就是利用BP的思想来求，可以参考：一文。
　　代价函数关于权值矩阵A的导数如下（拓扑和非拓扑时结果是一样的，因为此时这2种情况下代价函数关于A是没区别的）：
　　非拓扑结构下代价函数关于s的导数如下：
　　拓扑Sparse coding下代价函数关于s的导数为：
　　关于本程序的一点注释：&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
如果按照上面公式的和我们的理解，A是由学习到的基向量构成，S为每个样本在该基分解下的系数。在这里表示前提下，可以这样定义：A为n*k维，其中的每一列表示的是训练出来的基向量，S是k*m,其中的每一列是对应输入样本的sparse coding分解系数，当然了，此时的X是n*m了。即每一列表示的是一个样本数据。如果我们反过来表示（虽然这样理解不对，这里只是用不同的表示方法矩阵而已），即A表示输入数据X的分解系数（即编码值），而S是原始数据集训练出来的基的构成的，那么此时关于A,S,X三者的维度可以这样定义和解释：现假设有m个样本X，每个样本是个n维的向量，即X为m*n维的矩阵，需要用sparse coding学习k个特征，使得代价函数值最小，则其中的A是m*k维的，A中的第i行表示第i个样本分解后的系数值，S是k*n维的，S的第i行表示所学习到的第i个基。当然了，在本次实验和以后类似情况下我们还是用正确的版本，即第一种表示。
在matlab中，右除矩阵A和右乘inv(A)虽然在定义上式一样的，但是两者运行的结果有可能不同，右除的精度要高些。
注意拓扑结构下代价函数对s导数公式中的最后一项是点乘符号，也就是矩阵中对应元素的相乘，如果弄成了普通的矩阵乘法则就很难通过gradient checking了。
本程序训练样本IMAGE原图片尺寸512*512，共10张，从这10张大图片中提取2w张8*8的小patch图片，这些图片部分显示如下：
　　一些Matlab函数：
　　circshift:
　　该函数是将矩阵循环平移的函数，比如说B = circshift(A,shiftsize)是将矩阵A按照shiftsize的方式左右平移，一般hiftsize为一个多维的向量，第一个元素表示上下方向移动（更准确的说是在第一个维度上移动，这里只是考虑是2维矩阵的情况，后面的类似），如果为正表示向下移，第二个元素表示左右方向移动，如果向右表示向右移动。
　　rndperm：
　　该函数是随机产生一个行向量，比如randperm(n)产生一个n维的行向量，向量元素值为1~n，随机选取且不重复。而randperm(n,k)表示产生一个长为k的行向量，其元素也是在1到n之间，不能有重复。
　　questdlg：
　　button = questdlg('qstring','title','str1','str2','str3',default)，这是一个对话框，对话框中的内容用qstring表示，标题为title，然后后面3个分别为对应yes,no,cancel按钮，最后的参数default为默认的对应按钮。
　　实验结果：
　　交叉优化参数中，给定s优化A时，由于A有直接的解析解，所以不需要通过lbfgs等优化算法求得，通过令代价函数对A的导函数为0，可以得到解析解为：
　　注意单位矩阵前一定要有个系数（即样本个数），不然在程序中直接用该方法求得的A是通过不了验证。
　　此时学习到的非拓扑结果为：
　　上面的结果有点难看，采用的是16*16大小的patch，而非8*8的。　　
　　采用cg优化，256个16*16大小的patch，其结果如下：
　　如果将patch改为8*8,121个特征点，结果如下（这个比较像样）：
　　如果用lbfgs，256个16*16的，其结果如下（效果很差，说明优化方法对结果有影响）：
　　实验部分代码及注释：
　　sparseCodeingExercise.m:
%% CS294A/CS294W Sparse Coding Exercise
Instructions
------------
This file contains code that helps you get started on the
sparse coding exercise. In this exercise, you will need to modify
sparseCodingFeatureCost.m and sparseCodingWeightCost.m. You will also
need to modify this file, sparseCodingExercise.m slightly.
% Add the paths to your earlier exercises if necessary
% addpath /path/to/solution
%%======================================================================
%% STEP 0: Initialization
Here we initialize some parameters used for the exercise.
addpath minF
numPatches = 20000;
% number of patches
numFeatures = 256;
% number of features to learn
patchDim = 16;
% patch dimension
visibleSize = patchDim * patchD %单通道灰度图，64维，学习121个特征
% dimension of the grouping region (poolDim x poolDim) for topographic sparse coding
poolDim = 3;
% number of patches per batch
batchNumPatches = 2000; %分成10个batch
lambda = 5e-5;
% L1-regularisation parameter (on features)
epsilon = 1e-5; % L1-regularisation epsilon |x| ~ sqrt(x^2 + epsilon)
gamma = 1e-2;
% L2-regularisation parameter (on basis)
%%======================================================================
%% STEP 1: Sample patches
images = load('IMAGES.mat');
images = images.IMAGES;
patches = sampleIMAGES(images, patchDim, numPatches);
display_network(patches(:, 1:64));
%%======================================================================
%% STEP 3: Iterative optimization
Once you have implemented the cost functions, you can now optimize for
the objective iteratively. The code to do the iterative optimization
using mini-batching and good initialization of the features has already
been included for you.
However, you will still need to derive and fill in the analytic solution
for optimizing the weight matrix given the features.
Derive the solution and implement it in the code below, verify the
gradient as described in the instructions below, and then run the
iterative optimization.
% Initialize options for minFunc
options.Method = 'cg';
options.display = 'off';
options.verbose = 0;
% Initialize matrices
weightMatrix = rand(visibleSize, numFeatures);%64*121
featureMatrix = rand(numFeatures, batchNumPatches);%121*2000
% Initialize grouping matrix
assert(floor(sqrt(numFeatures)) ^2 == numFeatures, 'numFeatures should be a perfect square');
donutDim = floor(sqrt(numFeatures));
assert(donutDim * donutDim == numFeatures,'donutDim^2 must be equal to numFeatures');
groupMatrix = zeros(numFeatures, donutDim, donutDim);%121*11*11
groupNum = 1;
for row = 1:donutDim
for col = 1:donutDim
groupMatrix(groupNum, 1:poolDim, 1:poolDim) = 1;%poolDim=3
groupNum = groupNum + 1;
groupMatrix = circshift(groupMatrix, [0 0 -1]);
groupMatrix = circshift(groupMatrix, [0 -1, 0]);
groupMatrix = reshape(groupMatrix, numFeatures, numFeatures);%121*121
if isequal(questdlg('Initialize grouping matrix for topographic or non-topographic sparse coding?', 'Topographic/non-topographic?', 'Non-topographic', 'Topographic', 'Non-topographic'), 'Non-topographic')
groupMatrix = eye(numFeatures);%非拓扑结构时的groupMatrix矩阵
% Initial batch
indices = randperm(numPatches);%1*20000
indices = indices(1:batchNumPatches);%1*2000
batchPatches = patches(:, indices);
fprintf('%6s%12s%12s%12s%12s\n','Iter', 'fObj','fResidue','fSparsity','fWeight');
for iteration = 1:200
iteration = 1;
error = weightMatrix * featureMatrix - batchP
error = sum(error(:) .^ 2) / batchNumP
%说明重构误差需要考虑样本数
fResidue = error;
num_batches = size(batchPatches,2);
R = groupMatrix * (featureMatrix .^ 2);
R = sqrt(R + epsilon);
fSparsity = lambda * sum(R(:));
%稀疏项和权值惩罚项不需要考虑样本数
fWeight = gamma * sum(weightMatrix(:) .^ 2);
%上面的那些权值都是随机初始化的
%10.4f\n', iteration, fResidue+fSparsity+fWeight, fResidue, fSparsity, fWeight)
% Select a new batch
indices = randperm(numPatches);
indices = indices(1:batchNumPatches);
batchPatches = patches(:, indices);
% Reinitialize featureMatrix with respect to the new
% 对featureMatrix重新初始化，按照网页教程上介绍的方法进行
featureMatrix = weightMatrix' * batchP
normWM = sum(weightMatrix .^ 2)';
featureMatrix = bsxfun(@rdivide, featureMatrix, normWM);
% Optimize for feature matrix
options.maxIter = 20;
%给定权值初始值，优化特征值矩阵
[featureMatrix, cost] = minFunc( @(x) sparseCodingFeatureCost(weightMatrix, x, visibleSize, numFeatures, batchPatches, gamma, lambda, epsilon, groupMatrix), ...
featureMatrix(:), options);
featureMatrix = reshape(featureMatrix, numFeatures, batchNumPatches);
weightMatrix = (batchPatches*featureMatrix')/(gamma*num_batches*eye(size(featureMatrix,1))+featureMatrix*featureMatrix');
[cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures, batchPatches, gamma, lambda, epsilon, groupMatrix);
display_network(weightMatrix);
　　sparseCodingWeightCost.m:
function [cost, grad] = sparseCodingWeightCost(weightMatrix, featureMatrix, visibleSize, numFeatures,
patches, gamma, lambda, epsilon, groupMatrix)
%sparseCodingWeightCost - given the features in featureMatrix,
computes the cost and gradient with respect to
the weights, given in weightMatrix
% parameters
weightMatrix
- the weight matrix. weightMatrix(:, c) is the cth basis
featureMatrix - the feature matrix. featureMatrix(:, c) is the features
for the cth example
visibleSize
- number of pixels in the patches
numFeatures
- number of features
- weight decay parameter (on weightMatrix)
- L1 sparsity weight (on featureMatrix)
- L1 sparsity epsilon
groupMatrix
- the grouping matrix. groupMatrix(r, :) indicates the
features included in the rth group. groupMatrix(r, c)
is 1 if the cth feature is in the rth group and 0
otherwise.
if exist('groupMatrix', 'var')
assert(size(groupMatrix, 2) == numFeatures, 'groupMatrix has bad dimension');
groupMatrix = eye(numFeatures);%非拓扑的sparse coding中，相当于groupMatrix为单位对角矩阵
numExamples = size(patches, 2);%测试代码时为5
weightMatrix = reshape(weightMatrix, visibleSize, numFeatures);%其实传入进来的就是这些东西
featureMatrix = reshape(featureMatrix, numFeatures, numExamples);
% -------------------- YOUR CODE HERE --------------------
% Instructions:
Write code to compute the cost and gradient with respect to the
weights given in weightMatrix.
% -------------------- YOUR CODE HERE --------------------
%% 求目标的代价函数
delta = weightMatrix*featureMatrix-
fResidue = sum(sum(delta.^2))./numE%重构误差
fWeight = gamma*sum(sum(weightMatrix.^2));%防止基内元素值过大
sparsityMatrix = sqrt(groupMatrix*(featureMatrix.^2)+epsilon);
fSparsity = lambda*sum(sparsityMatrix(:)); %对特征系数性的惩罚值
cost = fResidue+fWeight+fS %目标的代价函数
cost = fResidue+fW
%% 求目标代价函数的偏导函数
grad = (2*weightMatrix*featureMatrix*featureMatrix'-2*patches*featureMatrix')./numExamples+2*gamma*weightM
grad = grad(:);
　　sparseCodingFeatureCost&.m:
function [cost, grad] = sparseCodingFeatureCost(weightMatrix, featureMatrix, visibleSize, numFeatures, patches, gamma, lambda, epsilon, groupMatrix)
%sparseCodingFeatureCost - given the weights in weightMatrix,
computes the cost and gradient with respect to
the features, given in featureMatrix
% parameters
weightMatrix
- the weight matrix. weightMatrix(:, c) is the cth basis
featureMatrix - the feature matrix. featureMatrix(:, c) is the features
for the cth example
visibleSize
- number of pixels in the patches
numFeatures
- number of features
- weight decay parameter (on weightMatrix)
- L1 sparsity weight (on featureMatrix)
- L1 sparsity epsilon
groupMatrix
- the grouping matrix. groupMatrix(r, :) indicates the
features included in the rth group. groupMatrix(r, c)
is 1 if the cth feature is in the rth group and 0
otherwise.
isTopo = 1;
L = size(groupMatrix,1);
[K M] = size(featureMatrix);
if exist('groupMatrix', 'var')
assert(size(groupMatrix, 2) == numFeatures, 'groupMatrix has bad dimension');
if(isequal(groupMatrix,eye(numFeatures)));
isTopo = 0;
groupMatrix = eye(numFeatures);
isTopo = 0;
numExamples = size(patches, 2);
weightMatrix = reshape(weightMatrix, visibleSize, numFeatures);
featureMatrix = reshape(featureMatrix, numFeatures, numExamples);
% -------------------- YOUR CODE HERE --------------------
% Instructions:
Write code to compute the cost and gradient with respect to the
features given in featureMatrix.
You may wish to write the non-topographic version, ignoring
the grouping matrix groupMatrix first, and extend the
non-topographic version to the topographic version later.
% -------------------- YOUR CODE HERE --------------------
%% 求目标的代价函数
delta = weightMatrix*featureMatrix-
fResidue = sum(sum(delta.^2))./numE%重构误差
fWeight = gamma*sum(sum(weightMatrix.^2));%防止基内元素值过大
sparsityMatrix = sqrt(groupMatrix*(featureMatrix.^2)+epsilon);
fSparsity = lambda*sum(sparsityMatrix(:)); %对特征系数性的惩罚值
cost = fResidue++fSparsity+fW%此时A为常量，可以不用
cost = fResidue++fS
%% 求目标代价函数的偏导函数
gradResidue = (-2*weightMatrix'*patches+2*weightMatrix'*weightMatrix*featureMatrix)./numE
% 非拓扑结构时：
if ~isTopo
gradSparsity = lambda*(featureMatrix./sparsityMatrix);
% 拓扑结构时
gradSparsity = lambda*groupMatrix'*(groupMatrix*(featureMatrix.^2)+epsilon).^(-0.5).*featureM%一定要小心最后一项是内积乘法
gradSparsity = lambda*groupMatrix'*(groupMatrix*(featureMatrix.^2)+epsilon).^(-0.5).*featureM%一定要小心最后一项是内积乘法
grad = gradResidue+gradS
grad = grad(:);
　　sampleIMAGES.m:
function patches = sampleIMAGES(images, patchsize,numpatches)
% sampleIMAGES
% Returns 10000 patches for training
% load IMAGES;
% load images from disk
%patchsize = 8;
% we'll use 8x8 patches
%numpatches = 10000;
% Initialize patches with zeros.
Your code will fill in this matrix--one
% column per patch, 10000 columns.
patches = zeros(patchsize*patchsize, numpatches);
%% ---------- YOUR CODE HERE --------------------------------------
Instructions: Fill in the variable called "patches" using data
from IMAGES.
IMAGES is a 3D array containing 10 images
For instance, IMAGES(:,:,6) is a 512x512 array containing the 6th image,
and you can type "imagesc(IMAGES(:,:,6))," to visualize
it. (The contrast on these images look a bit off because they have
been preprocessed using using "whitening."
See the lecture notes for
more details.) As a second example, IMAGES(21:30,21:30,1) is an image
patch corresponding to the pixels in the block (21,21) to (30,30) of
for imageNum = 1:10%在每张图片中随机选取1000个patch，共10000个patch
[rowNum colNum] = size(images(:,:,imageNum));
for patchNum = 1:2000%实现每张图片选取1000个patch
xPos = randi([1,rowNum-patchsize+1]);
yPos = randi([1, colNum-patchsize+1]);
patches(:,(imageNum-1)*2000+patchNum) = reshape(images(xPos:xPos+patchsize-1,yPos:yPos+patchsize-1,...
imageNum),patchsize*patchsize,1);
%% ---------------------------------------------------------------
% For the autoencoder to work well we need to normalize the data
% Specifically, since the output of the network is bounded between [0,1]
% (due to the sigmoid activation function), we have to make sure
% the range of pixel values is also bounded between [0,1]
% patches = normalizeData(patches);
%% ---------------------------------------------------------------
function patches = normalizeData(patches)
% Squash data to [0.1, 0.9] since we use sigmoid as the activation
% function in the output layer
% Remove DC (mean of images).
patches = bsxfun(@minus, patches, mean(patches));
% Truncate to +/-3 standard deviations and scale to -1 to 1
pstd = 3 * std(patches(:));
patches = max(min(patches, pstd), -pstd) /%因为根据3sigma法则，95%以上的数据都在该区域内
% 这里转换后将数据变到了-1到1之间
% Rescale from [-1,1] to [0.1,0.9]
patches = (patches + 1) * 0.4 + 0.1;
　　实验总结：
　　拓扑结构的Sparse coding未完成，跑出来没有效果，还望有人指导下。
　　已解决非拓扑下的Sparse coding，那时候出现的问题是因为在代价函数中，重构误差那一项没有除样本数（下面博文回复中网友给的提示），导致代价函数，导数，以及A的解析解都相应错了。
　　但是拓扑Sparse Coding依旧没有训练出来，因为训练过程中代价函数的值不是递减的，而是基本无规律。
　　基本解决了拓扑下的Sparse coding。以前训练不出特征来主要原因是在sampleIMAGES.m里没有将最后的patches归一化注释掉（个人猜测：采样前的大图片是经过白化了的，所以如果后面继续用那个带误差的归一化，可能引入更大的误差，导致给定的样本不适合Sparse coding），另外就是根据群里网友@地皮菜的提示，将优化算法由lbfgs改为cg就可以得出像样的结果。由此可见，不同的优化算法对最终的结果也是有影响的。
　　参考资料：
阅读(...) 评论()
阿萨德发斯蒂芬